Sommaire

Introduction
Analyses exploratoire et descriptive
Sélection des variables
Modèles de régression
Annexes

1. Introduction

La grande quantité d’informations créées et collectées chaque jour présente de nombreux avantages d’un point de vue statistique et analytique, mais peut parfois être un inconvénient dans le sens où toutes les données ne sont pas pertinentes dans l’étude d’un phénomène. Le développement des données massives (Big Data), a donc été accompagné par l’essor de nouvelles méthodes de sélection de variables qui visent à identifier les informations pertinentes de celles redondantes ou inutiles. En effet, un modèle est considéré comme explicite, interprétable lorsqu’il est parcimonieux et contient 4-5 variables qui pourront alors être interprétées librement.

La quantité de méthodes de sélection de variables est elle aussi importante, et l’objectif de ce dossier est d’en expérimenter certaines dans le but d’appliquer une modélisation par les Moindres Carrés Ordinaires (MCO), la plus parcimonieuse possible. Notre objectif est d’expliquer le cours du Baltic Dry Index (BDI) entre février 2007 et décembre 2018, avec des données mensuelles. Le BDI est un indice couvrant le coût du fret maritime ; il fournit pour cela une évaluation du prix à payer pour transporter par voie de mer différentes matières premières spécifiques telles que le charbon, les céréales ou les minerais.

2. Analyses exploratoire

Commençons dans une première partie par explorer nos données pour les comprendre puis appliquer certaines modifications pour les rendre exploitables pour la suite de l’analyse.

2.1 Importation

# Import des données
library(readxl)
df <- read_excel("./Données/database project fret.xlsx")
df <- data.frame(df)

Ainsi importée notre base contient 153 lignes et 29 colonnes dont une de temps avec les dates, une de la variable à expliquer et 27 variables explicatives. De manière à rendre la base exploitable pour l’analyse nous retirons les valeurs manquantes, et nous mettons le bon format à la variable de temps.

2.2 Analyse de la stationnarité

Nous analysons à présent la stationnarité des variables puisque celle-ci est nécessaire pour l’application des futures méthodes : elle implique que la série soit indépendante du temps, c’est à dire avec une même espérance mathématique et une même variance pour chaque observation. Nous utiliserons le test Augmented Dickey-Fuller (ADF) qui suit la règle de décision suivante :

\(H_0\) : existence d’une racine unitaire, la série n’est pas stationnaire
\(H_1\) : la série est stationnaire

# on met en série temporelle
library(tseries)
ts_fret <- ts(data = df_fret[,2:29], start=c(2007,01,01),frequency=12)

pvalue <- 0
for (i in 1:ncol(ts_fret)){
  pvalue[i] <- adf.test(ts_fret[,i])$p.value }
test_ADF <- as.data.frame(pvalue)

resultat <- 0   # Oui stationnaire
for (i in 1:nrow(test_ADF)){
  if (test_ADF[i,"pvalue"]<0.05)
    resultat[i] <- "Oui"
  else
    resultat[i] <- "Non"
}
test_ADF <- round(test_ADF,3)
test_ADF <- cbind(test_ADF,resultat)
rownames(test_ADF) <- paste(colnames(ts_fret))

# On affiche les résultats du test
test_ADF <- data.frame(test_ADF)
DT::datatable(test_ADF)

Par la sortie du test on voit que seules 6 variables sont considérées stationnaires sur 28. Cependant, comme visible sur le graphique ci-dessous où nous regardons l’évolution d’une variable prise aléatoirement, malgré le fait que le test ADF la considère stationnaire, nous voyons que ses fluctuations ne sont pas indépendantes du temps. Par précaution, nous decidons donc de stationnariser toutes les variables de notre base via une fonction lag qui les différencie toutes une fois (excepté la date) :

# différenciation en passant par le lag
df_fret[,2:29] <- df_fret[,2:29] %>% mutate(df_fret[,2:29] - lag(df_fret[,2:29]))
df_fret <- na.omit(df_fret)

2.3 Analyse des outliers

Nous continuons notre démarche de manipulations des données dans le but de les rendre exploitables, en identifiant et retirant les observations atypiques de la base. Certaines valeurs sont effectivement considérées comme atypiques, aberrantes par leur forte ou très faible valeur : il est important de les déceler car elles peuvent biaiser l’analyse. La moyenne et la médiane en sont de bons exemples ; une forte valeur va tirer la moyenne vers le haut tandis que la véritable valeur pour laquelle 50% de l’échantillon a plus et 50% a moins, sera plus faible car non influencée par cet outlier. En séries temporelles, ces derniers sont le plus souvent liés à des chocs exogènes tels que guerres, crises, ou encore changements politiques.

Selon la nature de la variable à laquelle nous souhaitons retirer les points atypiques la méthode diffère, nous distinguons donc 2 méthodes : une pour les indices financiers et l’autre pour les variables temporelles non financières.

a) Pour les indices

ts_fret <- ts(data = df_fret[,2:29], start=c(2007,01,01),frequency=12)
    # sur indice financier 
library(robustbase)
library(PerformanceAnalytics)
df_clean1 <- Return.clean(ts_fret[,c(1,2,4,5,7,11,15,16,26,27,28)], method = "boudt")

b) Pour les autres

library(tsoutliers)
df_clean2<-data.frame()

a <- df_fret[,-c(1,2,3,5,6,8,12,16,17,27,28,29)]
j<-0
for (n in names(a)){
  y<-ts(a[,c(as.character(n))])
  fit<-tso(y)
  i<-fit$yadj
  if(j==0){df_clean2<-i} else{df_clean2<-data.frame(df_clean2,i)}
  j<-j+1
}
names(df_clean2)<-names(a)

Après avoir rassemblé les variables débarrassées des outliers en un seul data frame, notre base est prête à l’emploi puisqu’elle est propre et les variables sont stationnarisées. Nous pouvons désormais nous pencher sur la visualisation des données pour comprendre grâce à des graphiques, des statistiques et des classifications, les liens entre les variables.

2.4 Graphiques

Dans un premier temps nous regardons les graphiques de certaines séries en niveau c’est à dire n’ayant subies aucune modification, de manière à comparer une variable considérée stationnaire par le test ADF et une considérée non stationnaire, puis nous regarderons l’effet de la stationnarisation sur l’évolution temporelle de ces deux variables.

On constate sur ce graphique à 4 arguments l’apparence des variables en niveau et lorsqu’elles ont subi une transformation. En rouge nous trouvons les variations en niveau d’une variable (“Crude_Oil_Prices_Brent”) classée comme stationnaire selon le test de la racine unitaire, et en bleu une variable en niveau aussi mais testée comme non stationnaire par le test ADF (“Capacity_Utilization”). Comme noté précédemment, on voit à travers ces graphiques que pour les variables considérées comme stationnaires par le test comme pour celles qui sont classées ‘non stationnaires’, on retrouve des tendances largement observables : à la fois sur une dépendance de la variance et/ou de la moyenne au temps. Il apparaît donc que le test statistique n’est pas fiable à 100% et il était plus prudent de stationnariser toutes les variables. Étant donné que la variable que nous expliquons est un cours boursier, les valeurs que nous avons obtenues en laguant celle-ci correspondent aux rentabilités de cet indice maritime. Aussi, si nous regardons les variables transformées en noir correspondant aux variables de la ligne supérieure, nous voyons que l’application d’un lag a bel et bien permis de les rendre stationnaires et donc exploitables. Le reste des graphiques sur les séries en niveau puis transformées est visible en annexes 1 et 2.

2.5 Statistiques descriptives

De nombreux packages sous R offrent des statistiques intéressantes et concrètes pour aider à comprendre les données. Nous confronterons les libraires sumarytools, DataExplorer et explore pour avoir l’analyse la plus complète possible.

	Baltic.Dry.Index	Business_Tendency_Surveys	Capacity_Utilization	Consumer_Sentiment	Container.Index	CPI	Crude.Oil.Prices..Brent.	Geopolitical.Risk	GEPU_current	GEPU_ppp	GISS_Temperature	Global.Economic.Policy	Gold_Princing_Index	Indice_Kilian	Industrial_Capacity	Industrial_Production	M2	News_Based_Policy_Uncert_Index	OECD_6NME_Industrial_Production	Petroleum_and_other_liquids_stocks_US	Spread	Taux_change_effectif_USD	Three_Component_Index	Trade_Weighted_USD	US_Ending_Stocks_Crude_Oil	VIX	world.trade..volume.	WTI
Mean	-7.1196589	0.0086862	0.0439876	0.0986014	0.3385330	0.3586420	1.4648215	0.1548356	0.2486528	1.0239145	-0.0006294	1.5946482	4.0498230	-0.6014748	0.0770021	0.6310961	48.5566693	0.8718842	0.2410656	1.215212e+03	0.0296775	0.0278403	0.5421065	0.0669000	6.894266e+02	-0.0257188	0.2806021	0.0170252
Std.Dev	656.4581705	0.2318855	0.3863986	3.6830213	0.6349491	0.5908415	4.9312831	32.3000748	23.1487096	25.5696148	0.1149294	29.0505524	45.2669116	22.5232269	0.1367976	2.0099995	29.7463368	41.4782801	0.6000363	1.581245e+04	0.1812535	1.2580521	24.7494860	1.3620267	1.178014e+04	3.8482404	1.1367261	6.2954028
Min	-2374.9606920	-0.9002421	-1.2425000	-10.0000000	-1.8448880	-1.9473910	-9.6600000	-87.2815308	-66.0240392	-74.2506137	-0.3700000	-112.0146010	-112.7548701	-72.1362351	-0.2961000	-4.1528815	-48.2250000	-135.4581375	-1.1952980	-3.512885e+04	-0.4192926	-3.1900000	-79.5132746	-3.3840623	-2.608395e+04	-10.2344737	-2.0590263	-20.0753991
Q1	-201.0000000	-0.1270776	-0.2010493	-2.2000000	0.1000000	0.0660000	-1.9795275	-15.1690370	-14.9257249	-15.6748095	-0.0700000	-15.7829616	-23.7750000	-14.2202300	0.0274000	-0.8340000	28.6000000	-21.1182327	-0.1853093	-9.208750e+03	-0.1010390	-0.8500000	-13.5552979	-0.8381975	-6.644000e+03	-2.2282900	-0.5547856	-2.9600000
Median	12.0000000	0.0178164	0.0378524	0.1000000	0.4000000	0.4490000	1.3600794	-1.4399523	0.6862719	2.8872247	0.0100000	1.3353041	3.4166667	0.5042700	0.1018000	0.5403000	47.9750000	-1.2262497	0.2609783	2.009200e+03	0.0224256	-0.0100000	-1.3853302	0.0874833	9.872500e+02	-0.2776190	0.1751966	0.9700000
Q3	274.0000000	0.1371510	0.3074945	2.8000000	0.7000000	0.7100000	5.0718012	14.1656227	13.3702425	14.6778127	0.0700000	13.9871300	36.0869565	15.3922590	0.1759000	1.9362000	63.9000000	21.7205963	0.5984157	1.371505e+04	0.1359873	0.8100000	14.6959763	0.9089964	8.262250e+03	1.2384897	1.0026139	3.8600000
Max	2084.0000000	0.6476482	1.1910050	8.1000000	1.8000000	2.2550000	13.7300000	103.4815413	70.1629117	78.6549286	0.3600000	112.0693264	150.5639312	69.4436060	0.3001000	6.1583000	162.4250000	137.8643419	1.9941843	3.613100e+04	0.5708677	3.3900000	80.5878525	5.0531346	3.239710e+04	14.6279644	3.7730078	14.2800000
MAD	355.8240000	0.2001973	0.3683555	3.4099800	0.4447800	0.4388496	5.2038180	21.8369623	19.9933128	21.9335881	0.1037820	22.0456390	46.4551370	21.9907416	0.1098607	2.0375372	26.5088880	32.3848538	0.6071236	1.729586e+04	0.1738166	1.2157320	20.8601815	1.3272837	1.122780e+04	2.5730780	1.1742136	5.3818380
IQR	466.0000000	0.2577544	0.4995219	4.8000000	0.6000000	0.6375000	6.8906633	29.0645407	27.0208029	30.2580642	0.1400000	29.3117050	58.7344720	29.2678000	0.1469000	2.7174212	34.9925000	42.5979977	0.7764162	2.278620e+04	0.2297041	1.5900000	27.3524017	1.7403302	1.465638e+04	3.4200565	1.5407147	6.7550000
CV	-92.2035989	26.6959113	8.7842714	37.3526274	1.8755897	1.6474411	3.3664738	208.6088918	93.0964993	24.9724126	-182.6099710	18.2175300	11.1775037	-37.4466699	1.7765439	3.1849340	0.6126107	47.5731519	2.4891003	1.301210e+01	6.1074474	45.1882036	45.6542867	20.3591292	1.708687e+01	-149.6275213	4.0510245	369.7700879
Skewness	-0.7121659	-0.4511541	0.0143366	-0.2310208	-1.0936282	-0.8735057	-0.0365857	0.3303020	-0.0734509	-0.0078405	-0.2312626	0.3657294	0.0426178	-0.2930919	-0.9785227	0.4405574	0.4606377	0.3593092	0.0927017	-2.914237e-01	0.3990852	0.1432087	0.3202576	0.3832184	9.493790e-02	1.1428916	0.3181931	-0.7820450
SE.Skewness	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	0.2027312	2.027312e-01	0.2027312	0.2027312	0.2027312	0.2027312	2.027312e-01	0.2027312	0.2027312	0.2027312
Kurtosis	3.7826961	2.6622277	0.5982618	0.0087747	3.0418242	3.5753337	-0.2825557	1.7626045	0.5917051	0.8203816	0.5533605	3.4018798	0.1957058	0.8520057	0.5785901	-0.0134365	1.4413528	1.7293429	0.0048461	-6.522712e-01	0.3201489	-0.1374408	1.2823305	0.5710535	-1.807981e-01	3.5261932	-0.2565568	0.9096755
N.Valid	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	143.0000000	1.430000e+02	143.0000000	143.0000000	143.0000000	143.0000000	1.430000e+02	143.0000000	143.0000000	143.0000000
Pct.Valid	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	100.0000000	1.000000e+02	100.0000000	100.0000000	100.0000000	100.0000000	1.000000e+02	100.0000000	100.0000000	100.0000000

Grâce à la libraire sumarytools nous constatons que les rentabilités du Baltic Dry Index s’étendent de -2374.96 à 2084 avec une moyenne de -7.12 pour une médiane de 12. En comparaison avec l’étendue de la série (i.e. l’écart entre les valeurs maximales et minimales), la différence entre la moyenne et la médiane est minime puisqu’il est de 19.12 et représente donc 0.43% de l’étendu global de la série - soit moins de 1%. Ainsi la suppression des valeurs atypiques a-t-elle été efficace pour lisser les valeurs et ne pas altérer les calculs.

La libraire DataExplorer nous offre elle un aperçu des observations manquantes sur la base initiale, pour laquelle on observe un taux de valeurs manquantes de 5.88% pour 21 d’entre elles, ce qui représente 9 observations. Nous avons donc procédé dès le début de l’analyse à la suppression des observations incomplètes, passant ainsi de 153 à 144 mensualités. Puis en ayant différencié nos variables nous avons perdu une observation de plus, notre base finale contient donc 143 dates pour 29 variables. De même, nous observons les histogrammes de distribution de deux variables ci-dessous, nous constatons alors que les distributions s’apparentent davantage à des lois normales lorsque les données sont nettoyées. Les données nettoyées correspondent aux variables pour lesquelles on a retiré les observations manquantes, que l’on a stationnarisées et qui ont été débarrassées des points atypiques. Le changement est clair pour ces 2 variables “Crude.Oil.Prices..Brent” et “M2” - les distributions de toutes les variables de la base avant et après transformation sont quant à elles disponibles en annexe n°3.

2.6 Méthodes de classification

Dans cette partie nous allons nous intéresser à la classification de nos données; cette phase d’analyse peut nous permettre d’avoir un premier aperçu sur les variables qui divisent au mieux le jeu, celles qui sont les plus pertinentes dans l’explication de la variable à expliquer. Ainsi nous commençons par construire un arbre sur le Baltic Dry Index, puis nous réaliserons une analyse de composantes principales (ACP) pour finir sur une méthode ayant pour objectif de mettre en avant des clusters, c’est à dire des groupements de nos variables.

Le fonctionnement de l’arbre est le suivant : à chaque noeud, l’algorithme considère les N variables et cherche celle qui divise le jeu de données de la manière la plus optimale possible en maximisant la diminution globale de l’impureté, c’est-à-dire en réduisant le plus possible l’erreur. Plus simplement, à chaque noeud l’algorithme considère les variables les plus influentes et divise à partir de celles-ci un ensemble en 2 groupes qui sont les plus hétérogènes entre eux, mais les plus homogènes en leur sein. Appliqué sur nos données, l’arbre nous permet d’avoir une première idée sur l’importance de chaque variable et des liens qui les relient entre elles.

On voit ainsi qu’en partant de la valeur moyenne du BDI, à savoir -7.12 et du jeu complet avec 143 observations, la variable qui dichotomise le jeu de données initial en 2 groupes les plus homogènes possible est celle de l’indice Kilian. Pour 88 observations de cette dernière, la valeur est supérieure à -5.9. Pour ces 88 dates, l’algorithme du CART les divise à nouveau en 2 groupes homogènes où l’un a une capacité industrielle inférieure à -0.036 et est composé de 12 observations, tandis que l’autre a une capacité industrielle supérieure ou égale à -0.036 et englobe 76 observations soit 86.36% de ce sous-échantillon. Pour les 55 observations pour lesquelles l’indice Kilian est inférieur à -5.9, c’est la variable “VIX” qui permet de diviser ce sous-échantillon en 2 groupes hétérogènes entre eux.

Nous constatons donc à chaque noeud une nouvelle division du jeu ; finalement grâce à cet arbre de décision nous nous rendons compte que les variables importantes sélectionnées sont celles de l’indice Kilian (importante à 29% dans l’explication d’Y), de l’indice VIX (à 11%), de la capacité industrielle à 8%, du taux de change effectif et des enquêtes de conjoncture à 7% chacune. Nous nous attendons à ce que les variables dichotomisant le jeu aux noeuds principaux et ayant une grande importance soient significatives dans l’explication de l’indice de fret maritime, et ainsi qu’elles soient sélectionnées par les méthodes que nous appliquerons plus tard pour construire des modèles parcimonieux en dernière partie.

Pour compléter l’analyse de classification des variables commencée avec l’arbre de décision, nous appliquons une analyse en composantes principales (ACP) qui nous permet de grouper les variables autour d’axes fictifs et ainsi de voir les liens qui existent entre elles. Sur le cercle de corrélations des variables, il est difficile d’identifier les variables ayant des caractéristiques communes puisque les noms sont difficilement lisibles. On remarque tout de même en groupement de variables le long de l’axe n°1 avec des cosinus carrés très proches selon le code couleur : il s’agit des variables ‘GEPU_current’, ‘Global.Economic.Policy’, ‘New_Based_Policy_Uncert_Index’, ‘GEPU_ppp’ et ‘Three_Component_Index’ qui semblent par leur nom axées globalement sur la politique économique. Par l’annexe n°4 où sont les pourcentages de variance expliquée par les axes fictifs, on voit que la répartition de la variance expliquée par les axes est plutôt égale ; la première et la deuxième dimensions expliquent ainsi entre 30 et 31% de la variance. En regardant les graphiques de contribution des variables aux dimensions 1 et 2, nous voyons que l’axe 1 regroupe principalement les variables sur la politique économique mondiale dont nous venons de noter la proximité, tandis que l’axe 2 est axé sur le commerce extérieur avec les taux de change, les prix aux échanges etc. L’ACP nous a permis de constater les regroupements de variables tels que :

Global.Economic.Policy, GEPU_current, GEPU_ppp, News_Based_Policy_Uncert_Index et Three_Component_Index
Trade_Weighted_USD, Gold_Princing_Index, Taux_change_effectif_USD et Crude.Oil.PRices..Brent

Dans un dernier temps nous cherchons à regrouper autrement les variables ; par des Clusters en utilisant le package ClustOfVar qui offre un regroupement des variables et non des observations contrairement aux autres packages sous R. Cet algorithme de clustering a pour objectif d’identifier des groupes d’observations ayant des caractéristiques similaires avec le même principe que l’arbre CART, c’est-à-dire où les observations dans un même groupe se ressemblent le plus possible mais où elles se démarquent dans les groupes différents.

Ainsi, on applique la fonction hclustvar() sur notre jeu de données, en regardant le dendogramme des variables : on observe alors quelques groupements intéressants tels que:

Global Economic Policy avec GEPU_ppp et GEPU_current dont les initiales signifient Global Economic Policy Uncertaintly (nous nous doutons de la redondance de ces variables)
World trade volume avec la production industrielle des pays de l’OCDE
Crude oil prices brent avec WTI qui est l’indice boursier du pétrole (le lien ici est très fort)
Trade Weighted USD avec Taux_change_effectif_USD qui, bien que n’ayant pas les mêmes valeurs, semblent être très proches

Sur le graphique de la stabilité des partitionnements on constate que la stabilité est assez constante bien qu’elle diminue légèrement selon les partitions jusqu’à 5 classes. Nous décidons de garder 3 clusters et de regarder les variables qui sont contenues dans chacune des classes.

Le premier groupe rassemble les variables “Crude.Oil.Prices..Brent”, “WTI”, “CPI” avec une forte corrélation positive avec la variable synthétique, et “Taux_change_effectif_USD” et “Trade_Weighted_USD” avec relation négative forte. Le cluster n°1 représente donc le commerce extérieur (et correspondrait ainsi à la dimension 2 de l’ACP). La deuxième classe rassemblant davantage des variables axées sur la production telles que “Container.Index”, “OECD_6NME_Industrial_Production”, “Capacity_Utilization”, correspondrait donc à la capacité productive du pays. Enfin, le cluster n°3 s’apparente à la dimension 1 de l’ACP puisqu’elle regroupe les variables de politique mondiale qui semblent extrêmement corrélées puisqu’elles ont des corrélations proches les unes les autres avec la variable fictive du cluster n°3, et qu’elles semblent mesurer les mêmes choses - du moins c’est ce qui paraît avec la proximité de leurs noms et de leurs coefficients.

L’analyse en clustering nous a ainsi permis d’alimenter le constat sur la redondance de certaines variables explicatives, que nous avions pu noter dès la visualisation de l’arbre de décision et de l’analyse en composantes principales. Nous comprenons grâce à cette analyse descriptive, l’importance de la sélection de variables avant modélisation pour éviter les problèmes de mutlicollinéarité.

2.7 Corrélation

DataExplorer::plot_correlation(df_fret_clean)

Aussi, après avoir regardé différentes classifications pour comprendre les liens entre nos variables, nous regardons les coefficients de corrélations (CC) que les lient entre elles. Nous voyons alors sur la matrice de corrélations que certaines variables ont un CC supérieur à 0.5, nous savons alors qu’elles ne pourront être intégrées dans un même modèle car leur dépendance peut être gênante pour les estimations. Aussi les principales corrélations sont les suivantes :

Global.Economic.Policy, GEPU_ppp et GEPU_current pour lesquelles nous avions déjà souligné un lien fort
News_Based_Policy_Uncert_index avec Three_Component_Index qui étaient regroupées dans les mêmes classes lors des visualisations de données
Crude.Oil.Prices..Brent et “WTI” qui comme souligné précédemment correspond à l’indice du pétrole ; la forte corrélation entre ces 2 variables n’est donc pas surprenante

Nous devinons alors aisément que les variables “GEPU_ppp”, “GEPU_current”, “News_Based_Policy_Uncert_index” et “WTI” seront certainement éliminées des explicatives en raison de leur corrélation forte voire de leur redondance avec d’autres variables. Finalement, avant de commencer la sélection des variables nous pouvons regarder ces quelques graphiques qui viennent confirmer tous nos constats et conclure sur la non-pertinence de certaines variables qui contiennent les mêmes informations que d’autres.

3. Sélection des variables

Si nous voulions modéliser le Baltic.Dry.Index nous pourrions le faire par un modèle de régression linéaire multivarié sur la base de la minimisation de la variance résiduelle pour l’estimation de ses coefficients. Le problème de telles modélisations réside dans l’instabilité des résultats de l’estimation, car aucune contrainte n’est appliquée aux paramètres du modèle. Autrement dit, plusieurs modèles peuvent aboutir à des solutions proches en termes de minimisation de l’erreur, mais de faibles changements dans les données peuvent produire des modèles très différents.

C’est pour cela que l’on a recours à diverses méthodes de régularisation afin de rétrécir l’espace de solution autour de 0 : c’est pourquoi nous procéderons à la normalisation des variables explicatives et au centrage de la variable à expliquer pour enlever la constante. Ces méthodes de rétrécissement de l’espace de solution pour empêcher les valeurs trop grandes des estimateurs \(\beta\) sont appelées “méthodes de shrinkage”.

Plusieurs méthodes existent alors pour restreindre l’espace de solutions :

Certaines méthodes ont une approche plutôt “économétrique”. Parmi les n variables explicatives, on va rechercher le plus petit sous-ensemble de variables qui explique notre Y.
- leaps and bounds de Furnival et Wilson, qui consiste à tester des modèles pour tous les arrangements possibles. Mais cela devient rapidement infaisable avec un grand nombre de variables.
- régression pas à pas (stepwise) c’est une méthode itérative qui peut se décliner en deux variantes, forward qui initialise un modèle constitué uniquement de la constante, puis ajoute séquentiellement variable par variable celles qui améliorent le mieux l’ajustement en termes de réduction de la variabilité résiduelle. Backward qui fonctionne sur le même principe que la sélection forward mais de manière inversée ; on initialise la procédure avec un modèle complet comprenant toutes les variables que l’on va retirer séquentiellement.

D’autres méthodes plus “statistiques”, utilisent une fonction de pénalisation avec un paramètre \(\lambda\) (fixé empiriquement), que l’on rajoute à la fonction coût (d’une régression) afin d’obtenir de meilleurs résultats lors de l’estimation par validation croisée. Un \(\lambda\) élevé pénalisera plus fortement la régression, augmentant le biais des estimateurs mais réduisant la variance. Tout l’enjeu consiste donc à trouver le juste \(\lambda\) qui réponde au mieux à ce compromis entre biais et variance. Autrement dit, il s’agit du compromis entre la complexité d’un modèle et sa capacité à se généraliser. Le modèle optimal est un modèle qui est capable de décrire le mieux les données d’apprentissage, et capable de faire les meilleures prévisions possible tout en restant robuste.

Avant de commencer le processus de séléction nous allons donc normaliser les variables explicatives et centrer notre variable à expliquer Baltic.Dry.Index, afin de retirer le problème de la constante lors de l’estimation.

y <- df_fret_clean %>%
  select(Baltic.Dry.Index) %>%
  scale(center = T, scale = F) %>%
  as.matrix()

x <- df_fret_clean %>%
  select(-c(1,2)) %>%      # on retire la variable à expliquer et la date
  as.matrix()

3.1 Approche économétrique : GETS

Théorie

La méthode GETS (general-to-specific), est une méthode de la famille des backward selections c’est-à-dire qui commence à partir d’un modèle général et qui, après estimation de plusieurs modèles, procède à l’élimination des variables non significatives seulement si les hypothèses sont respectées à chaque spécification. L’objectif final de cette méthode est de trouver un modèle qui n’omette pas de variables explicatives importantes (bonne spécification) et qui retire toutes les variables superflues (parcimonieux).

L’application de la méthode GETS sur notre jeu de données se fera avec le package gets. Cette méthode se déroule en plusieurs étapes :

La première consiste à estimer un modèle non restreint général, ce qui dans notre cas sera un modèle autorégressif d’ordre 1 AR(1)
La deuxième consiste à réaliser des tests sur les résidus afin de vérifier le respect des hypothèses :
- Tests d’autocorrélation sur les résidus standardisés (Ljung et Box, 1979)
- Tests d’hétéroscédasticité conditionnelle sur les résidus standardisés au carré (Ljung et Box, 1979)
- Test de normalité (Jarque et Bera, 1980)

Application

library(stargazer)
library(lgarch)
library(gets)

tbase <- df_fret_clean[,-1]
class(tbase[,2:28]) #tibble
## [1] "data.frame"

#class(newbase[,2:28])
mX = data.matrix(tbase[,2:28])
#str(mX)

# ARX model 

AR_BDI <- arx(tbase$Baltic.Dry.Index, mc = T, ar = 1, mxreg = mX[,c(1:27)], vcov.type = "white", normality.JarqueB=T)
AR_BDI
## 
## Date: Sat Dec 05 13:11:49 2020 
## Dependent var.: y 
## Method: Ordinary Least Squares (OLS) 
## Variance-Covariance: White (1980) 
## No. of observations (mean eq.): 142 
## Sample: 2 to 143 
## 
## Mean equation:
## 
##                                              coef   std.error  t-stat   p-value
## mconst                                 4.9622e+01  1.1449e+02  0.4334  0.665530
## ar1                                   -1.3766e-01  1.0267e-01 -1.3408  0.182678
## Container.Index                        1.3730e+02  9.2627e+01  1.4823  0.141032
## Global.Economic.Policy                 3.0502e+00  2.1088e+00  1.4464  0.150838
## Geopolitical.Risk                      3.6672e-01  1.2423e+00  0.2952  0.768392
## WTI                                    1.8186e+01  1.1167e+01  1.6285  0.106209
## CPI                                   -3.3201e+02  1.0524e+02 -3.1548  0.002058
## Gold_Princing_Index                    1.0907e+00  1.0845e+00  1.0058  0.316680
## Indice_Kilian                          1.8820e+01  2.1407e+00  8.7917 1.889e-14
## Three_Component_Index                 -7.1201e+00  5.9855e+00 -1.1896  0.236712
## News_Based_Policy_Uncert_Index         6.3595e-01  3.9137e+00  0.1625  0.871206
## VIX                                    7.6749e+00  1.3740e+01  0.5586  0.577542
## Crude.Oil.Prices..Brent.               5.2406e+00  1.3494e+01  0.3884  0.698479
## world.trade..volume.                  -3.3557e+01  3.8405e+01 -0.8738  0.384091
## Business_Tendency_Surveys              2.9765e+01  2.7419e+02  0.1086  0.913747
## Capacity_Utilization                  -3.1741e+01  1.2688e+02 -0.2502  0.802921
## Consumer_Sentiment                    -2.7053e+01  1.2092e+01 -2.2373  0.027231
## GISS_Temperature                       3.8689e+02  4.0305e+02  0.9599  0.339162
## GEPU_current                           4.1743e+00  5.4467e+00  0.7664  0.445049
## GEPU_ppp                              -9.6701e+00  4.7846e+00 -2.0211  0.045636
## Industrial_Production                  7.5867e+00  1.8144e+01  0.4181  0.676638
## Industrial_Capacity                    5.3654e+02  4.0119e+02  1.3374  0.183784
## M2                                    -1.5660e-01  1.4549e+00 -0.1076  0.914476
## OECD_6NME_Industrial_Production        2.9587e+00  8.8774e+01  0.0333  0.973472
## Petroleum_and_other_liquids_stocks_US  2.5050e-03  2.8320e-03  0.8845  0.378284
## Spread                                -8.0665e+01  3.4091e+02 -0.2366  0.813386
## Trade_Weighted_USD                    -1.9859e+02  7.1883e+01 -2.7627  0.006693
## Taux_change_effectif_USD               1.5159e+02  7.3288e+01  2.0685  0.040878
## US_Ending_Stocks_Crude_Oil             3.4892e-03  3.6557e-03  0.9545  0.341880
##                                          
## mconst                                   
## ar1                                      
## Container.Index                          
## Global.Economic.Policy                   
## Geopolitical.Risk                        
## WTI                                      
## CPI                                   ** 
## Gold_Princing_Index                      
## Indice_Kilian                         ***
## Three_Component_Index                    
## News_Based_Policy_Uncert_Index           
## VIX                                      
## Crude.Oil.Prices..Brent.                 
## world.trade..volume.                     
## Business_Tendency_Surveys                
## Capacity_Utilization                     
## Consumer_Sentiment                    *  
## GISS_Temperature                         
## GEPU_current                             
## GEPU_ppp                              *  
## Industrial_Production                    
## Industrial_Capacity                      
## M2                                       
## OECD_6NME_Industrial_Production          
## Petroleum_and_other_liquids_stocks_US    
## Spread                                   
## Trade_Weighted_USD                    ** 
## Taux_change_effectif_USD              *  
## US_Ending_Stocks_Crude_Oil               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostics and fit:
## 
##                    Chi-sq df    p-value
## Ljung-Box AR(2)    1.1175  2 5.7194e-01
## Ljung-Box ARCH(1) 21.4183  1 3.6923e-06
## Jarque-Bera       34.6182  2 3.0392e-08
##                             
## SE of regression   484.86136
## R-squared            0.56373
## Log-lik.(n=142)  -1065.09782

On peut constater que seules quelques variables sont significatives dans le General Unrestreicted Model GUM, et si on regarde les tests sur les résidus on voit qu’ils sont non corrélés au seuil de risque de 5% - par contre ils ne valident pas le test d’homoscédasticité ni celui de la normalité.

Dans cette deuxième étape nous estimons le modèle en appliquant l’algorithme GETS qui éliminera de manière successive (backward elimination) les variables non significatives du GUM à partir d’un seuil de significativité fixé à 5%. Comme la première étape, la validité de chaque suppression est vérifiée par rapport aux tests de diagnostic des résidus. L’algorithme prendra fin lorsqu’il aura sélectionné, parmi les modèles terminaux, la spécification la mieux ajustée selon le critère d’ajustement du BIC.

#getsm_BDI <- getsm(AR_BDI)
#getsm_BDI

L’estimation par la méthodologie Gets ne fonctionne pas si nous ne lui disons pas de ne pas appliquer l’hypothèse d’homoscédasticité des résidus. C’est l’une des limites que nous pouvons constater de cette méthode; il peut arriver que les tests de diagnostics (autocorrélation et hétéroscédasticité conditionnelle) soient trop restrictifs et ne sélectionnent aucune variable ou alors un très petit nombre. C’est pourquoi nous allons réestimer en relâchant la contrainte d’homoscédasticité des résidus (ARCH).

    # GETS modelling without ARCH test
getsm_BDI2 <- getsm(AR_BDI, arch.LjungB=NULL)
getsm_BDI2
## 
## Date: Sat Dec 05 13:11:51 2020 
## Dependent var.: y 
## Method: Ordinary Least Squares (OLS) 
## Variance-Covariance: White (1980) 
## No. of observations (mean eq.): 142 
## Sample: 2 to 143 
## 
## GUM mean equation:
## 
##                                       reg.no. keep        coef  std.error
## mconst                                      1    0  4.9622e+01 1.1449e+02
## ar1                                         2    0 -1.3766e-01 1.0267e-01
## Container.Index                             3    0  1.3730e+02 9.2627e+01
## Global.Economic.Policy                      4    0  3.0502e+00 2.1088e+00
## Geopolitical.Risk                           5    0  3.6672e-01 1.2423e+00
## WTI                                         6    0  1.8186e+01 1.1167e+01
## CPI                                         7    0 -3.3201e+02 1.0524e+02
## Gold_Princing_Index                         8    0  1.0907e+00 1.0845e+00
## Indice_Kilian                               9    0  1.8820e+01 2.1407e+00
## Three_Component_Index                      10    0 -7.1201e+00 5.9855e+00
## News_Based_Policy_Uncert_Index             11    0  6.3595e-01 3.9137e+00
## VIX                                        12    0  7.6749e+00 1.3740e+01
## Crude.Oil.Prices..Brent.                   13    0  5.2406e+00 1.3494e+01
## world.trade..volume.                       14    0 -3.3557e+01 3.8405e+01
## Business_Tendency_Surveys                  15    0  2.9765e+01 2.7419e+02
## Capacity_Utilization                       16    0 -3.1741e+01 1.2688e+02
## Consumer_Sentiment                         17    0 -2.7053e+01 1.2092e+01
## GISS_Temperature                           18    0  3.8689e+02 4.0305e+02
## GEPU_current                               19    0  4.1743e+00 5.4467e+00
## GEPU_ppp                                   20    0 -9.6701e+00 4.7846e+00
## Industrial_Production                      21    0  7.5867e+00 1.8144e+01
## Industrial_Capacity                        22    0  5.3654e+02 4.0119e+02
## M2                                         23    0 -1.5660e-01 1.4549e+00
## OECD_6NME_Industrial_Production            24    0  2.9587e+00 8.8774e+01
## Petroleum_and_other_liquids_stocks_US      25    0  2.5050e-03 2.8320e-03
## Spread                                     26    0 -8.0665e+01 3.4091e+02
## Trade_Weighted_USD                         27    0 -1.9859e+02 7.1883e+01
## Taux_change_effectif_USD                   28    0  1.5159e+02 7.3288e+01
## US_Ending_Stocks_Crude_Oil                 29    0  3.4892e-03 3.6557e-03
##                                          t-stat   p-value    
## mconst                                 0.433429  0.665530    
## ar1                                   -1.340791  0.182678    
## Container.Index                        1.482346  0.141032    
## Global.Economic.Policy                 1.446384  0.150838    
## Geopolitical.Risk                      0.295188  0.768392    
## WTI                                    1.628476  0.106209    
## CPI                                   -3.154802  0.002058 ** 
## Gold_Princing_Index                    1.005760  0.316680    
## Indice_Kilian                          8.791674 1.889e-14 ***
## Three_Component_Index                 -1.189562  0.236712    
## News_Based_Policy_Uncert_Index         0.162495  0.871206    
## VIX                                    0.558596  0.577542    
## Crude.Oil.Prices..Brent.               0.388362  0.698479    
## world.trade..volume.                  -0.873780  0.384091    
## Business_Tendency_Surveys              0.108556  0.913747    
## Capacity_Utilization                  -0.250156  0.802921    
## Consumer_Sentiment                    -2.237268  0.027231 *  
## GISS_Temperature                       0.959886  0.339162    
## GEPU_current                           0.766378  0.445049    
## GEPU_ppp                              -2.021071  0.045636 *  
## Industrial_Production                  0.418141  0.676638    
## Industrial_Capacity                    1.337382  0.183784    
## M2                                    -0.107635  0.914476    
## OECD_6NME_Industrial_Production        0.033328  0.973472    
## Petroleum_and_other_liquids_stocks_US  0.884538  0.378284    
## Spread                                -0.236613  0.813386    
## Trade_Weighted_USD                    -2.762731  0.006693 ** 
## Taux_change_effectif_USD               2.068469  0.040878 *  
## US_Ending_Stocks_Crude_Oil             0.954476  0.341880    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostics:
## 
##                    Chi-sq df    p-value
## Ljung-Box AR(2)    1.1175  2 5.7194e-01
## Ljung-Box ARCH(1) 21.4183  1 3.6923e-06
## Jarque-Bera       34.6182  2 3.0392e-08
## 
## Paths searched: 
## 
## path 1 : 1 15 11 24 26 16 23 5 12 21 13 19 14 18 25 2 8 6 28 3 -7 -22 
## path 2 : 2 23 24 11 15 26 16 21 13 5 12 1 19 14 18 25 8 6 28 3 -7 -22 
## path 3 : 3 16 13 23 11 24 15 5 26 12 25 21 14 18 8 22 2 6 19 28 -1 -7 
## path 4 : 4 24 23 26 15 16 21 13 5 12 11 1 19 14 18 25 2 8 6 28 -3 -20 -7 -22 
## path 5 : 5 24 15 11 23 16 26 21 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 6 : 6 24 26 1 23 16 11 15 12 5 21 14 25 8 18 2 13 19 28 3 -7 -22 
## path 7 : 8 23 15 24 16 11 21 26 13 5 12 1 14 25 18 19 2 6 28 3 -7 -22 
## path 8 : 10 23 15 24 26 16 5 21 13 12 1 19 14 29 2 18 28 8 6 22 -3 -7 
## path 9 : 11 24 15 23 26 16 5 21 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 10 : 12 15 24 23 11 26 16 5 13 21 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 11 : 13 24 11 15 23 26 16 21 5 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 12 : 14 15 16 23 26 11 24 1 5 21 12 13 19 18 25 2 8 6 28 3 -7 -22 
## path 13 : 15 24 23 11 26 16 5 21 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 14 : 16 24 15 23 11 26 5 21 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 15 : 18 24 23 26 15 13 11 16 21 5 12 1 19 14 25 2 8 6 28 3 -7 -22 
## path 16 : 19 23 24 15 11 16 5 21 13 26 12 1 14 18 25 2 8 6 28 3 -7 -22 
## path 17 : 21 24 11 23 15 16 26 5 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 18 : 22 24 15 11 26 13 23 5 12 16 25 3 8 14 21 18 2 6 19 28 -1 -7 
## path 19 : 23 24 15 11 26 16 5 21 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 20 : 24 15 23 11 26 16 5 21 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 21 : 25 24 23 11 15 26 13 21 16 5 12 1 14 18 19 2 8 6 28 3 -7 -22 
## path 22 : 26 24 15 23 11 16 5 21 13 12 1 19 14 18 25 2 8 6 28 3 -7 -22 
## path 23 : 29 24 11 15 23 5 16 26 12 13 21 1 19 14 18 2 28 8 6 22 -3 7 
## 
## Terminal models: 
## 
## spec 1 : 4 7 9 10 17 20 22 27 29 
## spec 2 : 1 4 7 9 10 17 20 27 29 
## spec 3 : 3 7 9 10 17 20 22 27 29 
## spec 4 : 3 4 7 9 11 17 20 25 27 
## spec 5 : 3 4 9 10 17 20 25 27 
## 
##          info(sc)      logl   n   k
## spec 1:    15.423 -1072.751 142   9
## spec 2:    15.427 -1072.991 142   9
## spec 3:    15.416 -1072.270 142   9
## spec 4:    15.435 -1073.604 142   9
## spec 5:    15.428 -1075.557 142   8
## 
## SPECIFIC mean equation:
## 
##                                   coef   std.error  t-stat   p-value    
## Container.Index             1.1458e+02  6.9103e+01  1.6581  0.099653 .  
## CPI                        -1.8857e+02  1.0232e+02 -1.8428  0.067579 .  
## Indice_Kilian               1.7850e+01  2.0050e+00  8.9029 3.523e-15 ***
## Three_Component_Index      -3.5914e+00  1.5236e+00 -2.3572  0.019871 *  
## Consumer_Sentiment         -2.8862e+01  1.1353e+01 -2.5422  0.012161 *  
## GEPU_ppp                   -3.1257e+00  1.7291e+00 -1.8077  0.072910 .  
## Industrial_Capacity         5.2904e+02  2.8158e+02  1.8788  0.062457 .  
## Trade_Weighted_USD         -1.1023e+02  3.5240e+01 -3.1279  0.002163 ** 
## US_Ending_Stocks_Crude_Oil  6.3810e-03  2.8836e-03  2.2129  0.028613 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostics and fit:
## 
##                    Chi-sq df    p-value
## Ljung-Box AR(2)    1.0336  2 5.9641e-01
## Ljung-Box ARCH(1) 20.1039  1 7.3346e-06
##                             
## SE of regression   475.29986
## R-squared            0.50657
## Log-lik.(n=142)  -1072.26959

Cette fois-ci la méthodologie GETS a fonctionné et nous pouvons observer à partir de la sortie précédente que le processus a réalisé 23 estimations au total. De ces 23, il en conserve 5 comme modèles terminaux et en sélectionne un seul. L’algorithme Gets a sélectionné un modèle contenant 9 variables explicatives, toutes significatives à 10%. Les variables explicatives sélectionnées expliquent environ 50% de la variance du Baltic.Dry.Index.

Malgré le fait d’avoir supprimé la contrainte d’homoscédasticité des résidus lors du processus de sélection, le modèle final ne remplit toujours pas cette hypothèse. Cela vient donc confirmer la limite de cette méthode de sélection de variables; c’est pourquoi d’autres méthodes moins contraignantes sont privilégiées.

Nous allons tout de même conserver les variables et leurs coefficients respectifs pour une comparaison générale entre méthodes de sélection.

    # GETS betas
gets_beta <- coef.arx(getsm_BDI2)
gets_beta
##            Container.Index                        CPI 
##               1.145794e+02              -1.885661e+02 
##              Indice_Kilian      Three_Component_Index 
##               1.785003e+01              -3.591435e+00 
##         Consumer_Sentiment                   GEPU_ppp 
##              -2.886183e+01              -3.125735e+00 
##        Industrial_Capacity         Trade_Weighted_USD 
##               5.290385e+02              -1.102281e+02 
## US_Ending_Stocks_Crude_Oil 
##               6.380973e-03

Fonction ‘isat’

Avant de passer à des méthodes de régressions pénalisées, le package gets propose la fonction isat() qui réalise la méthode de saturation d’indicateurs (IS) permettant la détection de ruptures et d’outliers.

Cette méthode permet 3 approches :
- Saturation d’indicateurs d’impulsion (IIS, impulse indicator saturation)
- Saturation d’indicateurs d’étape (SIS, step-indicator saturation)
- Saturation d’indicateurs de tendance (TIS, trend-indicator saturation)

Ci-dessous sont visibles les estimations par MCO via ces 3 approches ainsi que la représentation graphique du meilleur modèle. Celui-ci combine à la fois des indicateurs d’étape et de tendance avec un \(R^2\)=0.70069.

## 
## Date: Sat Dec 05 13:12:01 2020 
## Dependent var.: y 
## Method: Ordinary Least Squares (OLS)
## Variance-Covariance: Ordinary 
## No. of observations (mean eq.): 143 
## Sample: 1 to 143 
## 
## SPECIFIC mean equation:
## 
##            coef std.error   t-stat   p-value    
## mconst   496.71    141.33   3.5146 0.0006052 ***
## iis12  -2919.18    434.11  -6.7245 4.939e-10 ***
## iis20  -1509.83    403.30  -3.7437 0.0002704 ***
## iis21  -1820.86    399.79  -4.5545 1.187e-05 ***
## iis28   1782.10    380.01   4.6896 6.782e-06 ***
## sis8    1385.29    343.28   4.0354 9.215e-05 ***
## sis12   2723.03    508.29   5.3572 3.683e-07 ***
## sis17  -3972.05    389.29 -10.2034 < 2.2e-16 ***
## tis9   -1014.50    167.22  -6.0668 1.310e-08 ***
## tis13   1334.49    173.64   7.6852 3.184e-12 ***
## tis26  -1130.76    276.52  -4.0893 7.506e-05 ***
## tis27    811.71    239.80   3.3849 0.0009399 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostics and fit:
## 
##                   Chi-sq df  p-value
## Ljung-Box AR(1)   2.6378  1 0.104350
## Ljung-Box ARCH(1) 5.1189  1 0.023667
## 
##                     Stat. p-value
## Jiao-Pretis Prop. 4.32952 0.00001
## Jiao-Pretis Count 4.00000 0.00619
##                             
## SE of regression   373.91835
## R-squared            0.70069
## Log-lik.(n=143)  -1044.04557

3.2 Régressions pénalisées :

3.2.1 Ridge regression

Théorie

La régression Ridge a été initialement proposée pour limiter l’instabilité liée à des variables explicatives trop corrélées entre elles. La pénalisation ridge va diminuer la distance entre les solutions possibles sur la base de la distance euclidienne. Autrement dit elle ajoute des contraintes sur les estimateurs en réduisant leur amplitude, et par conséquent rendra faibles certains coefficients de l’estimation pour les variables moins pertinentes.

Quand le paramètre \(\lambda\) est proche de 0 on s’approche de la solution "classique", non pénalisée, et quand \(\lambda\) tend vers l’infini, la pénalisation est telle que tous les paramètres valent 0. Lorsque \(\lambda\) augmente, on augmente donc le biais de la solution mais on diminue sa variance. Afin d’optimiser le compromis biais-variance nous avons recours aux :

critères d’information (IC, information criteria)
critères de validation croisée (CV, cross-validation)

L’inconvénient de cette méthode est qu’elle ne sélectionne pas de variables explicatives car elles sont toutes présentent dans le modèle final ce qui nous empêche de conclure sur l’importance de chaque variable dans l’explication de la variable à expliquer.

Application

Nous allons à présent appliquer la régression ridge à nos données, en paramétrant à la fois la cross-validation à 10 estimations, et en définissant le nombre de \(\lambda\) qui sera testé (100), dont les valeurs seront sélectionnées aléatoirement parmi une plage de données renseignée allant de -3 à 5.

Le graphique ci-dessus nous présente les erreurs de prévisions par cross-validation pour toutes les valeurs de \(\lambda\) initialisées précédemment. La zone délimitée par les traits en pointillés correspond aux valeurs de \(\lambda\) pour lesquelles le Mean-Squared Error est minimisé.

L’étape suivante consiste à conserver les \(\lambda\) sélectionnés par la première estimation, et de ré-estimer le modèle à partir de ceux-ci : nous obtenons alors l’ensemble des coefficients (biaisés) de la régression ridge.

# Best lambda obtained from CV (lambda.min) - other possibility: lambda.1se
lambda_cv <- ridge_cv$lambda.min

# Evaluation of the final model with the selected lambda
model_cv <- glmnet(x, y, alpha = 0, lambda = lambda_cv, standardize = T)
model_cv
## 
## Call:  glmnet(x = x, y = y, alpha = 0, lambda = lambda_cv, standardize = T) 
## 
##   Df %Dev Lambda
## 1 27 49.8  215.4

## 27 x 1 sparse Matrix of class "dgCMatrix"
##                                                  s0
## Container.Index                        5.889561e+01
## Global.Economic.Policy                 4.599752e-01
## Geopolitical.Risk                      4.536502e-01
## WTI                                    1.031416e+01
## CPI                                   -1.118767e+02
## Gold_Princing_Index                    7.169629e-01
## Indice_Kilian                          1.221624e+01
## Three_Component_Index                 -2.081608e+00
## News_Based_Policy_Uncert_Index        -7.158984e-01
## VIX                                   -8.799664e-01
## Crude.Oil.Prices..Brent.              -8.986741e-01
## world.trade..volume.                  -2.247892e+01
## Business_Tendency_Surveys              1.777214e+02
## Capacity_Utilization                   5.156107e+01
## Consumer_Sentiment                    -2.243766e+01
## GISS_Temperature                       2.648191e+02
## GEPU_current                          -2.916771e-01
## GEPU_ppp                              -2.202959e+00
## Industrial_Production                 -2.640317e+00
## Industrial_Capacity                    3.810598e+02
## M2                                     1.377506e-01
## OECD_6NME_Industrial_Production       -2.079982e+01
## Petroleum_and_other_liquids_stocks_US  2.694537e-03
## Spread                                -1.365161e+02
## Trade_Weighted_USD                    -5.638689e+01
## Taux_change_effectif_USD               3.707335e+00
## US_Ending_Stocks_Crude_Oil             3.653594e-03

3.2.2 Lasso regression

Théorie

La méthode du LASSO (Least Absolute Shrinkage and Selection Operator) en comparaison à la régression ridge, va diminuer la distance des solutions possibles sur la base de la distance de Manhatan et non euclidienne. C’est une forme de pénalisation qui a la capacité d’éliminer les variables inutiles en fixant des coefficients à 0. Cela permet de simplifier les modèles en les rendant parcimonieux. Le LASSO ne fait donc pas de régression linéaire classique mais une régression régularisée qui rend nuls certains coefficients de l’estimation.

L’avantage du cette méthode est qu’elle réalise une vraie sélection de variables en tronquant à 0 les variables “inutiles”. En ce qui concerne ses limites en revanche, la régression Lasso pose problème en présence de variables corrélées car elle choisit aléatoirement une des deux variables et met l’autre à 0. De plus, avant toute estimation il est nécessaire de normaliser les données car cette méthode est sensible aux échelles des différentes variables explicatives dans son processus de sélection.

Application

De la même façon que pour ridge, nous allons lui indiquer les paramètres avant de lancer l’estimation. Puis nous conserverons les \(\lambda\) qui minimisent l’erreur et nous estimerons le modèle final à partir de ceux-ci.

# LASSO regression
# 10-folds CV to choose lambda
lambdas_to_try <- 10^seq(-3, 5, length.out = 100)


# alpha = 1, implementation of Lasso regression
lasso_cv <- cv.glmnet(x, y, alpha = 1, lambda = lambdas_to_try, standardize = T, nfolds = 10) # choix du meilleur lambda parmi 100

# Figures of lambdas
plot(lasso_cv)

Comme pour RIDGE, le graphique ci-dessus nous présente les erreurs de prévision par cross-validation pour les 100 valeurs de \(\lambda\) choisies aléatoirement.

# Best lambda obtained from CV (lambda.1se) - other possibility: lambda.min
lambda_cv <- lasso_cv$lambda.1se

# Evaluation of the final model with the selected lambda
model_cv <- glmnet(x, y, alpha = 1, lambda = lambda_cv, standardize = T)

## 27 x 1 sparse Matrix of class "dgCMatrix"
##                                             s0
## Container.Index                       .       
## Global.Economic.Policy                .       
## Geopolitical.Risk                     .       
## WTI                                   .       
## CPI                                   .       
## Gold_Princing_Index                   .       
## Indice_Kilian                         9.223744
## Three_Component_Index                 .       
## News_Based_Policy_Uncert_Index        .       
## VIX                                   .       
## Crude.Oil.Prices..Brent.              .       
## world.trade..volume.                  .       
## Business_Tendency_Surveys             .       
## Capacity_Utilization                  .       
## Consumer_Sentiment                    .       
## GISS_Temperature                      .       
## GEPU_current                          .       
## GEPU_ppp                              .       
## Industrial_Production                 .       
## Industrial_Capacity                   .       
## M2                                    .       
## OECD_6NME_Industrial_Production       .       
## Petroleum_and_other_liquids_stocks_US .       
## Spread                                .       
## Trade_Weighted_USD                    .       
## Taux_change_effectif_USD              .       
## US_Ending_Stocks_Crude_Oil            .

Après réestimation du modèle définitif à partir des paramètres choisis, on voit que la méthode ne conserve qu’une seule variable explicative pour le modèle et est donc très parcimonieuse, il s’agit de Indice_Kilian. Nous notons par ailleurs que cette variable est celle qui apparaissait la plus importante lors de la classification de notre base par l’arbre de décision, ce qui vient donc confirmer notre hypothèse selon laquelle les variables divisant le jeu de données aux noeuds principaux seraient certainement retenues pour la suite de l’analyse.

3.2.3 Weighted fusion Lasso regressions

Théorie

Les méthodes aLASSO et WF Fusion entre autres permettent de pallier au problème de sélection de variables via l’ajout de poids ; la régression Weighted Fusion (WF) rectifie le problème de corrélations entre les variables explicatives par l’ajout de données. La “fusion pondérée” peut potentiellement incorporer la redondance des informations entre les variables corrélées pour l’estimation et la sélection des variables. Elle est également utile lorsque le nombre de variables est supérieur au nombre d’observations. La WF Lasso permet ainsi d’améliorer la sélection des variables et la précision des prévisions.

Application

# Weighted fusion Lasso regressions

# Initialization of parameters
gamma=0.5
mu=0.1

# Compute pxQ matrix.
cor.mat <- cor(x)
abs.cor.mat <- abs(cor.mat)
sign.mat <- sign(cor.mat) - diag(2, nrow(cor.mat))
Wmat <- (abs.cor.mat^gamma - 1 * (abs.cor.mat == 1))/(1 -abs.cor.mat * (abs.cor.mat != 1))
weight.vec <- apply(Wmat, 1, sum)
fusion.penalty <- -sign.mat * (Wmat + diag(weight.vec))

# Compute Cholesky decomposition
R<-chol(fusion.penalty, pivot = TRUE)

# Transform Weighted Fusion in a Lasso issue
p<-dim(x)[2]
xstar<-rbind(x,sqrt(mu)*R)
ystar<-c(y,rep(0,p))

# Apply Lasso .
fit_wfusion<-glmnet(xstar,ystar)
fit_cv_wfusion<-cv.glmnet(xstar,ystar)

par(mfrow=c(1,2))
plot(fit_wfusion,xvar ="lambda",label = FALSE,col=T,xlab=expression(log(lambda)))
plot(fit_cv_wfusion,xlab=expression(log(lambda)))

min(fit_cv_wfusion$cvm)
## [1] 225837.8
lambda.opt<-fit_cv_wfusion$lambda.min

# We obtain the estimator of beta_weighted-fusion
wfusion_beta<-predict(fit_wfusion,type="coef",s=lambda.opt)
show(wfusion_beta)
## 28 x 1 sparse Matrix of class "dgCMatrix"
##                                                   1
## (Intercept)                            1.231102e+01
## Container.Index                        .           
## Global.Economic.Policy                 .           
## Geopolitical.Risk                      .           
## WTI                                    .           
## CPI                                    .           
## Gold_Princing_Index                    .           
## Indice_Kilian                          1.480575e+01
## Three_Component_Index                 -2.039285e+00
## News_Based_Policy_Uncert_Index         .           
## VIX                                    .           
## Crude.Oil.Prices..Brent.               .           
## world.trade..volume.                   .           
## Business_Tendency_Surveys              .           
## Capacity_Utilization                   .           
## Consumer_Sentiment                    -1.550831e+01
## GISS_Temperature                       .           
## GEPU_current                           .           
## GEPU_ppp                              -1.789590e+00
## Industrial_Production                  .           
## Industrial_Capacity                    .           
## M2                                     .           
## OECD_6NME_Industrial_Production        .           
## Petroleum_and_other_liquids_stocks_US  9.231731e-04
## Spread                                 .           
## Trade_Weighted_USD                    -5.833307e+01
## Taux_change_effectif_USD               .           
## US_Ending_Stocks_Crude_Oil             2.486618e-03

Avec cette méthode de séléction nous conservons 7 variables explicatives.

3.2.4 Elastic-Net regression

Théorie

En règle général, la régression ridge donne de meilleurs résultats que la LASSO lorsque les variables sont corrélées entre elles. Mais si on utilise ridge, on ne peut pas réduire le nombre de variables. Pour sortir de ce dilemme, il existe la régression Elastic-Net qui ajoute une pénalisation Ridge en plus de celle de LASSO et donc combine les deux approches.

Cette régression permet à la fois de gérer les problèmes liés aux corrélations entre variables explicatives (ridge), et de retirer les variables non pertinentes pour l’estimation d’un modèle final parcimonieux (LASSO).

Le paramètre \(\alpha\) est compris entre 0 et 1, il va permettre de définir l’équilibre entre ridge et LASSO. Lorsque \(\alpha = 1\) nous auront à faire à LASSO et pour \(\alpha = 0\) à la régression ridge. En effet, quand \(\alpha\) augmente de 0 à 1 avec \(\lambda\) fixe, le nombre de variables retirées du modèle augmente de façon monotone.

Application

Avant de procéder à la régression Elastic-Net nous devons donc approximer les valeurs des paramètres \(\alpha\) et \(\lambda\) qui minimiseraient l’erreur de cross-validation. Le code suivant permet cette recherche des hyperparamètres.

library(doParallel)
library(foreach)
# Elastic-Net regression
# Choose alpha sequencially with 0 < alpha < 1

lambdas_to_try <- 10^seq(-3, 5, length.out = 100)
a <- seq(0.1, 0.9, 0.05)
search <- foreach(i = a, .combine = rbind) %dopar% {
cv <- glmnet::cv.glmnet(x, y,  alpha = i, lambda = lambdas_to_try, standardize = T, nfolds = 10)
  data.frame(cvm = cv$cvm[cv$lambda == cv$lambda.1se], lambda.1se = cv$lambda.1se, alpha = i)
}

# Implementation of EN regression
elasticnet_cv <- search[search$cvm == min(search$cvm), ]
elasticnet_cv$alpha
## [1] 0.55
elasticnet_cv$lambda.1se
## [1] 178.865
elasticnet_cv$cvm
## [1] 295168
# Figures of lambdas

La valeur de \(\alpha\) est 0.55
La valeur de \(\lambda\) est 178.86
La moyenne des erreurs de cross-validation est 300931.7

Le graphique ci-dessus représente les erreurs de prévision par cross-validation pour les 100 valeurs de \(\lambda\) choisies aléatoirement entre -3 et 5.

Une fois ces valeurs connues, nous réalisons l’estimation en précisant ces paramètres et nous récupérons comme pour les autres estimations les coefficients \(\beta\) des variables explicatives conservées.

# Best lambda obtained from CV (lambda.1se) - other possibility: lambda.min
lambda_cv <- elasticnet_cv$lambda.1se

# Evaluation of the final model with the selected lambda
model_cv <- glmnet(x, y, lambda = elasticnet_cv$lambda.1se, alpha = elasticnet_cv$alpha)

## 27 x 1 sparse Matrix of class "dgCMatrix"
##                                                s0
## Container.Index                         .        
## Global.Economic.Policy                  .        
## Geopolitical.Risk                       .        
## WTI                                     .        
## CPI                                     .        
## Gold_Princing_Index                     .        
## Indice_Kilian                          11.2629417
## Three_Component_Index                  -0.4564160
## News_Based_Policy_Uncert_Index          .        
## VIX                                     .        
## Crude.Oil.Prices..Brent.                .        
## world.trade..volume.                    .        
## Business_Tendency_Surveys               .        
## Capacity_Utilization                    .        
## Consumer_Sentiment                     -3.8184658
## GISS_Temperature                        .        
## GEPU_current                            .        
## GEPU_ppp                               -0.5883527
## Industrial_Production                   .        
## Industrial_Capacity                     .        
## M2                                      .        
## OECD_6NME_Industrial_Production         .        
## Petroleum_and_other_liquids_stocks_US   .        
## Spread                                  .        
## Trade_Weighted_USD                    -31.3270532
## Taux_change_effectif_USD                .        
## US_Ending_Stocks_Crude_Oil              .

La régression Elastic-Net montre alors 5 variables importantes à garder pour les estimations, qui correspondent aux indices 7, 8, 15, 18 et 25 des variables explicatives.

3.2.5 Bridge regression

Théorie

En ce qui concerne la régression Bridge, elle vient comme la régression Elastic-Net faire un mixte entre les régressions Ridge et LASSO mais avec cette fois un paramètre \(q\) qui vient agir sur \(\lambda\) et contrôler le degré de biais affecté lors de l’estimation des coefficients. Lorsque \(q = 2\) nous avons à faire à une régression Ridge et lorsque \(q = 1\) à une régression LASSO. Pour avoir un modèle parcimonieux il est préférable de paramétrer \(q\) entre 0 et 1. C’est pouquoi nous avons fixé q à 0.5 dans notre cas pratique.

Application

# Bridge regression
# 10-folds CV to choose lambda
# seq(-3, 5, length.out = 100) : random sequence of 100 numbers betwwene -3 and 5
lambdas_to_try <- 10^seq(-3, 5, length.out = 100)

# choix du meilleur lambda parmi 100
bridge_cv <- cv.bridge(x, y, q=0.5, lambda = lambdas_to_try, nfolds = 10)
# Figures of lambdas
plot(bridge_cv)

# Best lambda obtained from CV (lambda.min) - other possibility: lambda.1se
lambda_cv <- bridge_cv$lambda.min

# Evaluation of the final model with the selected lambda
model_cv <- bridge(x, y, q=0.5, lambda = lambda_cv)
summary(model_cv)
##        Length Class  Mode   
## betas  28     -none- numeric
## lambda  1     -none- numeric

##                                       68926.1210434971
## (Intercept)                               1.017544e+01
## Container.Index                           0.000000e+00
## Global.Economic.Policy                    0.000000e+00
## Geopolitical.Risk                         0.000000e+00
## WTI                                       4.753113e-09
## CPI                                       0.000000e+00
## Gold_Princing_Index                       0.000000e+00
## Indice_Kilian                             1.691749e+01
## Three_Component_Index                     0.000000e+00
## News_Based_Policy_Uncert_Index            0.000000e+00
## VIX                                       0.000000e+00
## Crude.Oil.Prices..Brent.                  0.000000e+00
## world.trade..volume.                      0.000000e+00
## Business_Tendency_Surveys                 0.000000e+00
## Capacity_Utilization                      0.000000e+00
## Consumer_Sentiment                        0.000000e+00
## GISS_Temperature                          0.000000e+00
## GEPU_current                              0.000000e+00
## GEPU_ppp                                  0.000000e+00
## Industrial_Production                     0.000000e+00
## Industrial_Capacity                       0.000000e+00
## M2                                        0.000000e+00
## OECD_6NME_Industrial_Production           0.000000e+00
## Petroleum_and_other_liquids_stocks_US     0.000000e+00
## Spread                                    0.000000e+00
## Trade_Weighted_USD                       -3.346070e-07
## Taux_change_effectif_USD                  0.000000e+00
## US_Ending_Stocks_Crude_Oil                0.000000e+00

Après récupération du \(\lambda\) optimal, la régression Bridge a sélectionné 3 variables explicatives hors constante avec des coefficients non nuls. Nous passons maintenant à d’autres techniques de sélection de variables.

3.2.6 SCAD

Théorie

Le problème des régressions présentées précédemment à savoir LASSO, Ridge, Elastic Net et Bridge est que lorsque le nombre de variables pertinentes est très faible par rapport au nombre d’observations, les variables peuvent être parasites pour la régression et biaisées. Une des solutions à ce problème souligné par Fan, Li et Zhang consiste à augmenter le paramètre de lissage \(\lambda\). Une autre réside dans de nouvelles approches telles que la régression SCAD (“smoothly clipped absolute deviations”), qui introduit une fonction de pénalité bloquant la valeur d’un coefficient associé à une variable pertinente. De cette manière, la pénalité SCAD conserve le taux de pénalisation (et le biais) du LASSO pour les petits coefficients, mais l’assouplit continuellement à mesure que la valeur absolue du coefficient augmente - sélectionnant alors des modèles plus parcimonieux avec un biais réduit. Le fonctionnement de la pénalité SCAD est le suivant :

la pénalité est linéaire en \(\beta\) quand \(|\beta|≤ \lambda\)
la pénalité est quadratique quand \(\lambda ≤ |\beta| ≤ a\lambda\)
la pénalité est constante quand \(|\beta| > a\lambda\)

Nous comprenons alors que la régression SCAD pénalise plus sévèrement que le LASSO les petits coefficients, en revanche elle pénalise moins les grands coefficients, et pallie ainsi au problème soulevé antérieurement.

Application

library(ncvreg)
fit_SCAD=ncvreg(x, y, penalty = c("SCAD"))
plot(fit_SCAD)
summary(fit_SCAD, lambda=30) #On choisi 30 comme valeur de Lambda.
## SCAD-penalized linear regression with n=143, p=27
## At lambda=30.0000:
## -------------------------------------------------
##   Nonzero coefficients         :  13
##   Expected nonzero coefficients:   7.26
##   Average mfdr (13 features)   :   0.559
## 
##                                         Estimate       z      mfdr Selected
## Indice_Kilian                          1.761e+01  9.7893   < 1e-04        *
## Trade_Weighted_USD                    -1.099e+02 -3.6945 0.0099352        *
## Three_Component_Index                 -4.465e+00 -2.7439 0.0755948        *
## CPI                                   -1.781e+02 -2.6754 0.0900837        *
## Consumer_Sentiment                    -2.835e+01 -2.6411 0.0994021        *
## Gold_Princing_Index                    8.192e-03  0.7050 0.7371999        *
## Petroleum_and_other_liquids_stocks_US  1.835e-04  0.8150 0.7408786        *
## Container.Index                        1.338e+01  0.9530 0.7644979        *
## Industrial_Capacity                    1.371e+02  1.2059 0.8624980        *
## WTI                                    3.246e+00  1.2481 0.8845475        *
## world.trade..volume.                  -5.988e+00 -0.9111 1.0000000        *
## GEPU_ppp                              -1.556e+00 -1.6369 1.0000000        *
## US_Ending_Stocks_Crude_Oil             4.177e-03  1.7831 1.0000000        *

# Validation croisé pour le meilleur lambda
cvfit_SCAD=cv.ncvreg(x, y, penalty = c("SCAD"))
plot(cvfit_SCAD)

# On attribue le meilleur lambda
lambda_SCAD <- cvfit_SCAD$lambda.min

#Modele finale
SCAD_Final=ncvreg(x, y, lambda=lambda_SCAD, alpha = 1)

scad_beta <- SCAD_Final$beta
scad_beta
##                                             47.5738
## (Intercept)                            5.891723e+01
## Container.Index                        0.000000e+00
## Global.Economic.Policy                 0.000000e+00
## Geopolitical.Risk                      0.000000e+00
## WTI                                    0.000000e+00
## CPI                                   -1.063199e+02
## Gold_Princing_Index                    0.000000e+00
## Indice_Kilian                          1.741565e+01
## Three_Component_Index                 -3.802980e+00
## News_Based_Policy_Uncert_Index         0.000000e+00
## VIX                                    0.000000e+00
## Crude.Oil.Prices..Brent.               0.000000e+00
## world.trade..volume.                   0.000000e+00
## Business_Tendency_Surveys              0.000000e+00
## Capacity_Utilization                   0.000000e+00
## Consumer_Sentiment                    -2.293098e+01
## GISS_Temperature                       0.000000e+00
## GEPU_current                           0.000000e+00
## GEPU_ppp                              -1.571102e+00
## Industrial_Production                  0.000000e+00
## Industrial_Capacity                    0.000000e+00
## M2                                     0.000000e+00
## OECD_6NME_Industrial_Production        0.000000e+00
## Petroleum_and_other_liquids_stocks_US  0.000000e+00
## Spread                                 0.000000e+00
## Trade_Weighted_USD                    -1.072109e+02
## Taux_change_effectif_USD               0.000000e+00
## US_Ending_Stocks_Crude_Oil             4.050266e-03

La SCAD nous donne donc 7 variables à conserver pour l’estimation finale.

Complément

Bien que nous ne l’appliquions pas, nous trouvions intéressant d’expliquer rapidement le fonctionnement de la régression «minimax concave penalty» (MCP) puisqu’il est très proche de celui de la méthode de SCAD. Comme elle, le MCP commence par appliquer le même taux de pénalisation que le LASSO et le relâche lentement à mesure que la valeur des coefficients augmente. Comparé au SCAD cependant, le MCP relâche la pénalité immédiatement tandis que pour le SCAD le taux continue d’augmenter avant de décroître. Autrement dit, la fonction de pénalité du MCP est la même que celle du SCAD mais est composée seulement de 2 paliers (sur les 3 présentés dans la méthode SCAD).

3.2.7 Adaptive Lasso regression

Théorie

Par rapport à la régression LASSO, la méthode Adaptive LASSO (“aLASSO”) résout le problème de sélection de variables. En effet, bien que la régression LASSO ait de nombreux avantages, ses estimateurs sont biaisés et ce biais ne disparaît pas nécessairement lorsque \({n \to \infty}\). Ce biais dans l’estimation du LASSO est d’environ \(\lambda\) pour les coefficients de régression élevés. Or, étant donné que le biais de l’estimation est déterminé par \(\lambda\), une approche pour réduire le biais du lasso consiste à utiliser l’approche de pénalité pondérée. C’est le principe de la régression aLASSO : les poids sont adaptés aux variables et sont tels que les variables avec de coefficients élevés ont des poids plus faibles - ce qui permet de réduire le biais d’estimation du LASSO tout en conservant sa propriété de parcimonie, voire même améliorer la précision de sélection des variables. Ainsi la régression ‘aLASSO’ résout le problème de sélection de variables du LASSO par l’ajout d’un paramètre de pondération des variables.

Application

# Adaptive Lasso regression

#LASSO
model_cv <- glmnet(x, y, alpha = 1, lambda = lambda_cv, standardize = T)
coef_lasso<-predict(model_cv,type="coef",s=lambda_cv)

# Weighted with gamma = 0.5
gamma=0.5
w0<-1/(abs(coef_lasso) + (1/length(y)))
poids.lasso<-w0^(gamma)

# Adaptive LASSO
fit_adalasso <- glmnet(x, y, penalty.factor =poids.lasso)
fit_cv_adalasso<-cv.glmnet(x, y,penalty.factor=poids.lasso)

# Figure of lambdas
plot(fit_cv_adalasso)

# Best lambda obtained from CV (lambda.1se) - other possibility: lambda.min
lambda_cv <- fit_cv_adalasso$lambda.1se

# Evaluation of the final model with the selected lambda
model_cv <- glmnet(x, y, alpha = 1, lambda = lambda_cv, standardize = T)

# Lasso betas
alasso_beta <- model_cv$beta
alasso_beta
## 27 x 1 sparse Matrix of class "dgCMatrix"
##                                               s0
## Container.Index                        .        
## Global.Economic.Policy                 .        
## Geopolitical.Risk                      .        
## WTI                                    .        
## CPI                                    .        
## Gold_Princing_Index                    .        
## Indice_Kilian                         10.4093068
## Three_Component_Index                  .        
## News_Based_Policy_Uncert_Index         .        
## VIX                                    .        
## Crude.Oil.Prices..Brent.               .        
## world.trade..volume.                   .        
## Business_Tendency_Surveys              .        
## Capacity_Utilization                   .        
## Consumer_Sentiment                     .        
## GISS_Temperature                       .        
## GEPU_current                           .        
## GEPU_ppp                               .        
## Industrial_Production                  .        
## Industrial_Capacity                    .        
## M2                                     .        
## OECD_6NME_Industrial_Production        .        
## Petroleum_and_other_liquids_stocks_US  .        
## Spread                                 .        
## Trade_Weighted_USD                    -0.3384348
## Taux_change_effectif_USD               .        
## US_Ending_Stocks_Crude_Oil             .

Nous pouvons constater que la méthode aLASSO par rapport à la méthode LASSO sélectionne une variable supplémentaire “Trade_Weighted_USD”, l“indice_kilian” restant dans l’estimation. Par ailleurs on voit que son coefficient est plus fort que lors de l’estimation LASSO, à savoir 10.41 contre 9.22.

3.2.8 Adaptive Elastic-Net

Théorie

La régression “Adaptive Elastic Net” (aEN) a été présentée dès 2009 par Zou et Zhang qui ont réfléchi au problème de colinéarité pour la sélection et l’estimation de modèles, lorsque le nombre de paramètres diverge de la taille de l’échantillon. La méthode Elastic-Net adaptif a donc été développée : elle combine les forces de la régularisation quadratique à celles du lasso pondéré de manière adaptative. L’aEN traite donc mieux le problème de colinéarité que les autres méthodes de sélection de variables, bénéficiant ainsi de performances d’échantillon fini améliorées.

Application

Pour cette méthode comme pour adaptive SCAD, MS-aEN et MS-aSCAD nous utiliserons les fonctions contenues dans le package “msaenet” sous R.


## adaptive EN
# Choose alpha sequencially with 0 < alpha < 1
a <- seq(0.1, 0.9, 0.05)
aenet.fit <- aenet(x, y, family = "gaussian", init = "enet", alphas = a, tune = "cv", nfolds = 10, seed = 1001)
print(aenet.fit)
## Call: aenet(x = x, y = y, family = "gaussian", init = "enet", alphas = a, 
##     tune = "cv", nfolds = 10, seed = 1001) 
##   Df      %Dev       Lambda
## 1  5 0.4632652 1.960817e+17
coef(aenet.fit)
##  [1]   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
##  [7]  16.807955  -3.514441   0.000000   0.000000   0.000000   0.000000
## [13]   0.000000   0.000000 -28.235398   0.000000   0.000000  -3.136224
## [19]   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
## [25] -79.752388   0.000000   0.000000

# Adaptive EN beta
aen_beta <- aenet.fit$beta
aen_beta
## 27 x 1 sparse Matrix of class "dgCMatrix"
##                                               s0
## Container.Index                         .       
## Global.Economic.Policy                  .       
## Geopolitical.Risk                       .       
## WTI                                     .       
## CPI                                     .       
## Gold_Princing_Index                     .       
## Indice_Kilian                          16.807955
## Three_Component_Index                  -3.514441
## News_Based_Policy_Uncert_Index          .       
## VIX                                     .       
## Crude.Oil.Prices..Brent.                .       
## world.trade..volume.                    .       
## Business_Tendency_Surveys               .       
## Capacity_Utilization                    .       
## Consumer_Sentiment                    -28.235398
## GISS_Temperature                        .       
## GEPU_current                            .       
## GEPU_ppp                               -3.136224
## Industrial_Production                   .       
## Industrial_Capacity                     .       
## M2                                      .       
## OECD_6NME_Industrial_Production         .       
## Petroleum_and_other_liquids_stocks_US   .       
## Spread                                  .       
## Trade_Weighted_USD                    -79.752388
## Taux_change_effectif_USD                .       
## US_Ending_Stocks_Crude_Oil              .
# Get the name of relevant variables
which(! coef(aenet.fit) == 0, arr.ind = TRUE)   
## [1]  7  8 15 18 25

msaenet.nzv(aenet.fit)      # Get Indices of Non-Zero Variables
## [1]  7  8 15 18 25
msaenet.fp(aenet.fit, 1:5)  # Get the Number of False Positive Selections
## [1] 5
msaenet.tp(aenet.fit, 1:5)  # Get the Number of True Positive Selections
## [1] 0

La méthode ‘Elastic-Net adaptée’ appliquée à notre jeu de données donne les mêmes résultats que l’EN basique : ces régressions sélectionnent toutes deux 5 variables qui sont “indice_kilian”, “Three_Component_Index”, “Consumer_Sentiment”, “GEPU_ppp” et “Trade_Weighted_USD”. En revanche, le coefficient estimé pour chacune des variables sélectionnées est différent, ils semblent ainsi plus proches de 0 dans le cas de la méthode EN ; le \(\beta\) positif de l’indice Kilian est inférieur (12.34 pour EN contre 16.81 pour aEN) mais les betas négatifs sont supérieurs (-5.14 pour EN contre -28.24 pour aEN par exemple pour Consumer_Sentiment).

Les effets de chaque variable sur l’indice maritime semblent ainsi renforcés dans l’aEN par rapport à la régression Elastic-Net classique.

3.2.9 Adaptive SCAD-Net

## Adaptive SCAD-Net
a <- seq(0.1, 0.9, 0.05)
asnet.fit <- asnet(x, y, family = "gaussian", init = "snet", alphas = a, tune = "cv", nfolds = 10, seed = 1001)
print(asnet.fit)
## Call: asnet(x = x, y = y, family = "gaussian", init = "snet", alphas = a, 
##     tune = "cv", nfolds = 10, seed = 1001) 
##   Df  Lambda Gamma Alpha
## 1 16 35.5994   3.7   0.9
coef(asnet.fit)
##  [1]   55.12391928    0.00000000    0.00000000   10.76063649 -206.56964559
##  [6]    0.00000000   12.36775674   -1.37704558   -0.07084923    0.00000000
## [11]    0.00000000  -30.78956391  311.77391548   36.24916050  -25.74926298
## [16]  368.90591491    0.00000000   -1.35859272    0.00000000  507.79130476
## [21]    0.00000000  -26.23396803    0.00000000 -187.28929018  -79.45735460
## [26]    0.00000000    0.00000000

# Adaptive SCAD beta
ascad_beta <- asnet.fit$beta
ascad_beta 
## 27 x 1 sparse Matrix of class "dgCMatrix"
##                    
##  [1,]   55.12391928
##  [2,]    .         
##  [3,]    .         
##  [4,]   10.76063649
##  [5,] -206.56964559
##  [6,]    .         
##  [7,]   12.36775674
##  [8,]   -1.37704558
##  [9,]   -0.07084923
## [10,]    .         
## [11,]    .         
## [12,]  -30.78956391
## [13,]  311.77391548
## [14,]   36.24916050
## [15,]  -25.74926298
## [16,]  368.90591491
## [17,]    .         
## [18,]   -1.35859272
## [19,]    .         
## [20,]  507.79130476
## [21,]    .         
## [22,]  -26.23396803
## [23,]    .         
## [24,] -187.28929018
## [25,]  -79.45735460
## [26,]    .         
## [27,]    .

# Get the name of relevant variables
which(! coef(asnet.fit) == 0, arr.ind = TRUE)   
##  [1]  1  4  5  7  8  9 12 13 14 15 16 18 20 22 24 25

msaenet.nzv(asnet.fit)      # Get Indices of Non-Zero Variables
##  [1]  1  4  5  7  8  9 12 13 14 15 16 18 20 22 24 25
msaenet.fp(asnet.fit, 1:5)  # Get the Number of False Positive Selections
## [1] 13
msaenet.tp(asnet.fit, 1:5)  # Get the Number of True Positive Selections
## [1] 3

Par rapport à toutes les autres méthodes de sélection de variables que nous avons appliquées jusqu’alors, la régression “Adaptive SCAD-Net” est la moins parcimonieuse puisqu’elle indique que 16 variables sont à garder dans les estimations. Jusqu’à présent, c’était la méthode automatique GETS qui sélectionnait le plus grand nombre de variables (9). Parmi les variables écartées par l’aSCAD-N on retrouve “Three_Component_Index” et “Global.Economic.Policy” dont nous avions constaté précédemment la forte corrélation avec les autres régresseurs.

3.2.10 Multi-Step Adaptive Elastic-Net

Théorie

Le “Multi-step Adaptive Elastic Net” (MSA-Enet) a été développé en 2015 par Nan Xiao et Qingsong Xu. Cette méthode vise à réduire le nombre de faux-positifs, tout en maintenant la précision de l’estimation. Cela est rendu possible en analysant les variables éliminées à chaque étape, pour mieux comprendre la structure des groupes de variables corrélées. Son fonctionnement est le suivant :

les poids sont initialisés en début de procédure
la fonction de poids de l’aEN est appliquée
les poids adaptifs sont mis à jour

L’algorithme itère jusqu’à trouver les bonnes variables à retenir, et ainsi le MSA-Enet pallie aux nombreux problèmes de sélection et de régression de variables.

Application

## Multi-Step Adaptive Elastic-Net
a <- seq(0.1, 0.9, 0.05)
msaenet.fit <- msaenet(x, y, family = "gaussian", init = "enet", alphas = a, tune = "cv", nfolds = 10, seed = 1001)
print(msaenet.fit)
## Call: msaenet(x = x, y = y, family = "gaussian", init = "enet", alphas = a, 
##     tune = "cv", nfolds = 10, seed = 1001) 
##   Df      %Dev       Lambda
## 1  5 0.4627897 7.142639e+17
coef(msaenet.fit)
##  [1]   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
##  [7]  16.566865  -3.388452   0.000000   0.000000   0.000000   0.000000
## [13]   0.000000   0.000000 -27.997838   0.000000   0.000000  -2.956046
## [19]   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
## [25] -80.573044   0.000000   0.000000

# MS Adaptive EN beta
msaen_beta <- msaenet.fit$beta
msaen_beta
## 27 x 1 sparse Matrix of class "dgCMatrix"
##                                               s0
## Container.Index                         .       
## Global.Economic.Policy                  .       
## Geopolitical.Risk                       .       
## WTI                                     .       
## CPI                                     .       
## Gold_Princing_Index                     .       
## Indice_Kilian                          16.566865
## Three_Component_Index                  -3.388452
## News_Based_Policy_Uncert_Index          .       
## VIX                                     .       
## Crude.Oil.Prices..Brent.                .       
## world.trade..volume.                    .       
## Business_Tendency_Surveys               .       
## Capacity_Utilization                    .       
## Consumer_Sentiment                    -27.997838
## GISS_Temperature                        .       
## GEPU_current                            .       
## GEPU_ppp                               -2.956046
## Industrial_Production                   .       
## Industrial_Capacity                     .       
## M2                                      .       
## OECD_6NME_Industrial_Production         .       
## Petroleum_and_other_liquids_stocks_US   .       
## Spread                                  .       
## Trade_Weighted_USD                    -80.573044
## Taux_change_effectif_USD                .       
## US_Ending_Stocks_Crude_Oil              .

# Get the name of relevant variables
which(! coef(msaenet.fit) == 0, arr.ind = TRUE) 
## [1]  7  8 15 18 25

msaenet.nzv(msaenet.fit)        # Get Indices of Non-Zero Variables
## [1]  7  8 15 18 25
msaenet.fp(msaenet.fit, 1:5)    # Get the Number of False Positive Selections
## [1] 5
msaenet.tp(msaenet.fit, 1:5)    # Get the Number of True Positive Selections
## [1] 0

De la même manière que pour les régressions Elastic_net et Adaptive Elastic-net, la méthode retient les 5 mêmes variables dont les coefficients sont très proches de ceux estimés par aEN avec un écart de moins d’un point pour chaque \(\beta\).

3.2.11 Multi-Step Adaptive SCAD-Net

## Multi-Step Adaptive SCAD-Net (Est-ce la même chose que MD aSCAD ?)
a <- seq(0.1, 0.9, 0.05)
msasnet.fit <- msasnet(x, y, family = "gaussian", init = "snet", alphas = a, tune = "cv", nfolds = 10, seed = 1001)
print(msasnet.fit)
## Call: msasnet(x = x, y = y, family = "gaussian", init = "snet", alphas = a, 
##     tune = "cv", nfolds = 10, seed = 1001) 
##   Df   Lambda Gamma Alpha
## 1 15 62.33687   3.7   0.9
coef(msasnet.fit)
##  [1]   45.7583105    0.0000000    0.0000000    9.4460446 -184.2326825
##  [6]    0.0000000   10.2391214   -0.5089388    0.0000000    0.0000000
## [11]    0.0000000  -30.2715829  408.5037788   40.8010678  -24.2412242
## [16]  331.2411405    0.0000000   -0.3806617    0.0000000  506.7569016
## [21]    0.0000000  -32.5289653    0.0000000 -185.0829700  -80.1143216
## [26]    0.0000000    0.0000000

# MS Adaptive SCAD beta
msascad_beta <- msasnet.fit$beta
msascad_beta
## 27 x 1 sparse Matrix of class "dgCMatrix"
##                   
##  [1,]   45.7583105
##  [2,]    .        
##  [3,]    .        
##  [4,]    9.4460446
##  [5,] -184.2326825
##  [6,]    .        
##  [7,]   10.2391214
##  [8,]   -0.5089388
##  [9,]    .        
## [10,]    .        
## [11,]    .        
## [12,]  -30.2715829
## [13,]  408.5037788
## [14,]   40.8010678
## [15,]  -24.2412242
## [16,]  331.2411405
## [17,]    .        
## [18,]   -0.3806617
## [19,]    .        
## [20,]  506.7569016
## [21,]    .        
## [22,]  -32.5289653
## [23,]    .        
## [24,] -185.0829700
## [25,]  -80.1143216
## [26,]    .        
## [27,]    .

# Get the name of relevant variables
which(! coef(msasnet.fit) == 0, arr.ind = TRUE)
##  [1]  1  4  5  7  8 12 13 14 15 16 18 20 22 24 25

msaenet.nzv(msasnet.fit)        # Get Indices of Non-Zero Variables
##  [1]  1  4  5  7  8 12 13 14 15 16 18 20 22 24 25
msaenet.fp(msasnet.fit, 1:5)    # Get the Number of False Positive Selections
## [1] 12
msaenet.tp(msasnet.fit, 1:5)    # Get the Number of True Positive Selections
## [1] 3

La méthode Adaptive SCAD appliquée précédemment s’est révélée être la moins parcimonieuse de toutes puisqu’elle a sélectionné 16 régresseurs. La même méthode avec plusieurs étapes (multistep) donne les mêmes résultats à peu de choses près : 15 variables sont sélectionnées et cette fois l’indice “WTI” est écarté - variable pour laquelle nous avions noté précédemment une forte corrélation aux autres facteurs explicatifs.

3.3 Autre approche :

Théorie

Nous avons ainsi vu qu’il existe un nombre important de méthodes de sélection de variables, chacune avec ses avantages et se complétant les unes les autres. De nombreuses autres méthodes existent encore pour répondre à la problématique d’un grand nombre de données, telles que les méthodes proximales, la PLS, la régression en composantes principales, l’algorithme FISTA, le ScoTLASS, le smooth-LASSO, le kernelLASSO, les algorithmes génétiques, les forêts aléatoires, les réseaux de neurones MONNA, la méthode BoLASSO etc. Dans une dernière sous-section, c’est à cette dernière méthode que nous nous intéresserons ; nous commencerons par expliquer son fonctionnement théorique (comme pour les autres méthodes de sélection de variables), avant de l’appliquer sur notre base de données.

Développée par Francis R.Bach en 2008, la méthode BoLASSO correspond à une régression LASSO avec une perturbation bootstrap. Le bootstrap est une méthode de rééchantillonage visant à évaluer la performance générale d’un modèle en divisant la base afin d’estimer les performances sur un ensemble de validation. Un échantillon Boostrap étant construit par un tirage avec remise de n-observations parmi n-observations, il correspond à un cas particulier de l’échantillon d’origine car il garde les mêmes proportions mais est distribué de manière aléatoire. L’estimateur Bolasso est construit à partir de M échantillons bootstraps indépendants et est défini comme l’intersection de ces M ensembles amenant ainsi a un estimateur consistant en termes de sélection.

La méthode de sélection de variables “BoLasso” conduit à une sélection cohérente du modèle sans la condition de consistance requise par une régression LASSO basique. Francis Bach a effectivement montré dans un document de recherche, que lorsque la décroissance spécifique du paramètre de régularisation est proportionnelle à \(n^−1/2\), où n est le nombre d’observations, alors le Lasso sélectionne toutes les variables qui devraient entrer dans le modèle (les variables pertinentes) avec une probabilité augmentant exponentiellement à n - alors qu’il sélectionne toutes les autres variables avec une probabilité strictement positive. Le problème avec cette propriété est que si plusieurs ensembles de données générés à partir de la même distribution étaient disponibles, alors toutes les variables pertinentes seraient toujours sélectionnées pour tous les ensembles de données, mais les variables non pertinentes entreraient dans les modèles de manière aléatoire.

L’objectif de la régression Lasso Bootstrap est alors de casser ce processus de probabilités de sélection non adapté à toutes les variables selon leur pertinence, en rééchantillonnant à partir du même ensemble de données unique. Le Lasso est alors exécuté pour plusieurs réplications bootstrap d’un échantillon donné, et l’intersection des supports des estimations faites sur ces échantillons bootstrappés conduit à une sélection cohérente du modèle.

Application

Son application n’a pas été évidente puisque de nombreuses fonctions existent pour traiter des problèmes avec BoLASSO mais les packages ne sont plus disponibles pour cette version de R (package “mht” par exemple), d’autres s’avèrent être pour des problèmes catégoriels (package “SparseLearner”), d’autres encore sont développés par des programmateurs et impliquent d’autres fonctions pas toujours accessibles (package “dmolitor/bolasso” ou “david26694/bolasso”). Finalement, c’est la libraire ‘mhtboot’ que nous avons appelée, nous avons utilisé la fonction mht.1sample() qui implémente des tests d’hypothèses multiples basés sur la distribution bootstrap des p-values. Nous fixons le nombre de bootstrap à 100 et regardons les variables alors sélectionnées.

#install.packages("mhtboot")
library(mhtboot)
bolasso <- mht.1sample(x, B=100, ncpus=1)
# on regarde les variables sélectionnées
bolasso$signal
## [1]  2  4 11 14 19 21 26 27

La régression LASSO bootstrappée retient 8 variables comme la méthode “SCAD”. Par ailleurs les variables sélectionnées sont toutes différentes de cette dernière puisqu’il s’agit de Gold_Princing_Index, News_Based_Policy_Uncert_Index, Capacity_Utilization, GISS_Temperature, Industrial_Capacity, OECD_6NME_Industrial_Production, US_Ending_Stocks_Crude_Oil et VIX.

3.4 Comparaison

variables	gets_beta	Ridge	LASSO	Elastic_Net	WF_LASSO	Bridge	SCAD	AdapLASSO	AdapEN	AdapSCAD	MS_AdapEN	MS_AdapSCAD
Container.Index	114	58.8956109933933	0	0	0	0	0	0	0	55	0	45
Global.Economic.Policy	0	0.459975236940332	0	0	0	0	0	0	0	0	0	0
Geopolitical.Risk	0	0.453650150154125	0	0	0	0	0	0	0	0	0	0
WTI	0	10.3141575015133	0	0	0	0	0	0	0	10	0	9
CPI	-189	-111.876683110123	0	0	0	0	-107	0	0	-207	0	-185
Gold_Princing_Index	0	0.716962857834375	0	0	0	0	0	0	0	0	0	0
Indice_Kilian	17	12.2162408077372	9	11	14	16	17	10	16	12	16	10
Three_Component_Index	-4	-2.08160751013487	0	-1	-3	0	-4	0	-4	-2	-4	-1
News_Based_Policy_Uncert_Index	0	-0.715898423954272	0	0	0	0	0	0	0	-1	0	0
VIX	0	-0.879966380653333	0	0	0	0	0	0	0	0	0	0
Crude.Oil.Prices..Brent.	0	-0.898674071294302	0	0	0	0	0	0	0	0	0	0
world.trade..volume.	0	-22.4789165226824	0	0	0	0	0	0	0	-31	0	-31
Business_Tendency_Surveys	0	177.721398777374	0	0	0	0	0	0	0	311	0	408
Capacity_Utilization	0	51.5610716447208	0	0	0	0	0	0	0	36	0	40
Consumer_Sentiment	-29	-22.4376628202976	0	-4	-16	0	-23	0	-29	-26	-28	-25
GISS_Temperature	0	264.819138692403	0	0	0	0	0	0	0	368	0	331
GEPU_current	0	-0.291677128739515	0	0	0	0	0	0	0	0	0	0
GEPU_ppp	-4	-2.20295913472772	0	-1	-2	0	-2	0	-4	-2	-3	-1
Industrial_Production	0	-2.64031736182336	0	0	0	0	0	0	0	0	0	0
Industrial_Capacity	529	381.059847983298	0	0	0	0	0	0	0	507	0	506
M2	0	0.137750608583767	0	0	0	0	0	0	0	0	0	0
OECD_6NME_Industrial_Production	0	-20.7998237348071	0	0	0	0	0	0	0	-27	0	-33
Petroleum_and_other_liquids_stocks_US	0	0.00269453742525823	0	0	0	0	0	0	0	0	0	0
Spread	0	-136.516070457929	0	0	0	0	0	0	0	-188	0	-186
Trade_Weighted_USD	-111	-56.3868918112195	0	-32	-59	-1	-108	-1	-80	-80	-81	-81
Taux_change_effectif_USD	0	3.7073345623538	0	0	0	0	0	0	0	0	0	0
US_Ending_Stocks_Crude_Oil	0	0.00365359416525889	0	0	0	0	0	0	0	0	0	0
Total	9	toutes	1	5	7	3	7	2	5	16	5	15

La table ci-dessus nous informe des variables sélectionnées par chaque méthode ainsi que des coefficients estimés pour ces dernières, pour la méthode GETS le coefficient de la variable “US_Ending_Stocks_Crude_Oil” étant de 0.006 il a été arrondi à 0 mais nous la conserverons pour les estimations. Il en va de même pour les méthodes WF_LASSO, Bridge et enfin SCAD. On remarque :

“Global.Economic.Policy”, “Geopolitical.Risk”, “Gold_Princing_Index”, “VIX”, “Crude. Oil. Prices..Brent”, “GEPU_current”, “Industrial_Production”, “M2”, “Petroleum_and_other_liquids_stocks_US”, “Taux_change_effectif_USD” et “US_Ending_Stocks_Crude_Oil” ne sont sélectionnées par aucune méthode appliquée dans notre étude (hors méthode Ridge évidemment). Elles représentent tout de même 40.74% des variables explicatives.
“WTI”, “News_Based_Policy_Uncert_Index”, “world_trade..volume.”, “Business_Tendency_Surveys”, “Capacity_Utilization”, “GISS_Temperature”, “OECD_6NME_Industrial_Production” et “Spread” sont sélectionnées seulement par 1 ou 2 méthodes parmi les 11 appliquées (toujours hors Rige). Cela représente 29.63 % de la base.
A contrario, les variables les plus souvent sélectionnées (au minimum 8 fois sur 11) sont les suivantes : “Indice_Kilian”, “Three_Component_Index”, “Consumer_Sentiment”, “Trade_Weighted_USD” et “GEPU_ppp”. La variable qui arrive “en tête de liste” est l’indice Kilian qui est sélectionné par les 11 méthodes différentes et qui, nous le rappelons, était la variable la plus importante d’après la régression CART.

Nos méthodes de classification avaient donc vu juste quant aux groupements, aux fortes corrélations entre les variables. Ainsi parmi les 5 variables que nous avions pu voir en groupe sur l’ACP (Global.Economic.Policy, GEPU_current, GEPU_ppp, News_Based_Policy_Uncert_Index et Three_Component_Index) on voit que seules 2 d’entre elles ont été retenues quand les autres ont été écartées par les méthodes de sélection ; on voit cela puisque les variables “Global.Economic.Policy”, “GEPU_current” et “News_Based_Policy_Uncert_Index” n’ont pas ou très peu (moins de 3 fois) été retenues dans les différentes régressions. De la même manière, nous avions pu noter de fortes relations entre les variables suivantes : Trade_Weighted_USD, Gold_Princing_Index, Taux_change_effectif_USD et Crude.Oil.Prices..Brent. On voit cependant que seule cette première a été gardée dans les méthodes de sélection.

Parmi les variables les plus souvent sélectionnées, on peut constater que la corrélation entre les variables “GEPU_ppp” et “Three_Component_Index” est proche de 0.5, mais mise à part cette corrélation, les variables retenues par la plupart des méthodes ne présentent pas de patterns de dépendance avec d’autres variables. Nous pouvons ainsi, dans une dernière partie, estimer par régression linéaire différents modèles : un pour chaque méthode avec les variables sélectionnées qui seront mises comme régresseurs de Y - c’est-à-dire le Baltic Dry Index.

4. Modèle de régression

Nous voyons finalement toutes les estimations faites à partir des variables sélectionnées par les méthodes économétriques (GETS) et pénalisées (toutes les autres).

Chacun de ces modèles a été estimé sur la base propre composée de 142 observations et non 143, car en ajoutant le BDI lagué comme variable explicative nous avons perdu une observation supplémentaire. La première chose que l’on peut regarder pour chacun des modèles est la qualité d’ajustement ajustée (\(R^2\)) ainsi que les critères d’informations AIC et BIC. Ces mesures nous permettent de les comparer et de choisir ainsi les plus pertinents ; tel que le \(R^2\) est maximisé et les critères sont minimisés. Le critère d’Akaike (AIC) mesure la qualité d’un modèle statistique tandis que le critère bayésien (BIC) en est dérivé et prend en compte la taille de l’échantillon (dans notre cas tous les modèles sont construits sur la même base donc nous prêterons surtout intérêt aux critères AIC et \(R^2\)). De même, nous considérons le \(R^2\) ajusté et non le basique puisque celui-ci augmente à mesure que des variables explicatives sont intégrées aux modèles ; or il arrive que certaines n’améliorent pas réellement la qualité du modèle.

Code des modèles

# Estimation des modèles :-

  #Modèle GETS:
var_gets <- df_beta %>% select(1,2) %>% filter(gets_beta != 0) %>% select(1)

mod_gets <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Container.Index + CPI + Indice_Kilian + Consumer_Sentiment + GEPU_ppp + Industrial_Capacity + Trade_Weighted_USD + Three_Component_Index + US_Ending_Stocks_Crude_Oil, data = df_fret_clean)

  #Modèle LASSO:
var_lasso <- df_beta %>% select(1,4) %>% filter(LASSO != 0) %>% select(1)

mod_lasso <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Indice_Kilian, data = df_fret_clean)

  #Modèle Elastic-Net:
var_en <- df_beta %>% select(1,5) %>% filter(Elastic_Net!= 0) %>% select(1)

mod_en <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Indice_Kilian + Three_Component_Index + Consumer_Sentiment + GEPU_ppp + Trade_Weighted_USD, data = df_fret_clean)


 #Modèle WF_LASSO

var_wfusion <- df_beta %>% select(1,6) %>% filter(WF_LASSO != 0) %>% select(1)

mod_wfusion <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Indice_Kilian + Three_Component_Index + Consumer_Sentiment + GEPU_ppp + Petroleum_and_other_liquids_stocks_US + Trade_Weighted_USD + US_Ending_Stocks_Crude_Oil, data = df_fret_clean)


 #Modèle Bridge

var_bridge <- df_beta %>% select(1,7) %>% filter(Bridge != 0) %>% select(1)

mod_bridge <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + WTI + Indice_Kilian + Trade_Weighted_USD, data = df_fret_clean)


 #Modèle SCAD

var_scad <- df_beta %>% select(1,8) %>% filter(SCAD != 0) %>% select(1)

mod_scad <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + CPI + Indice_Kilian + Three_Component_Index + Consumer_Sentiment + GEPU_ppp + Trade_Weighted_USD + US_Ending_Stocks_Crude_Oil, data = df_fret_clean)


  #Modèle AdapLASSO:
var_alasso <- df_beta %>% select(1,9) %>% filter(AdapLASSO != 0) %>% select(1)

mod_alasso <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Indice_Kilian + Trade_Weighted_USD, data = df_fret_clean)


  #Modèle AdapEN:
var_aen <- df_beta %>% select(1,10) %>% filter(AdapEN != 0) %>% select(1)

mod_aen <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Indice_Kilian + Three_Component_Index + Consumer_Sentiment + GEPU_ppp + Trade_Weighted_USD, data = df_fret_clean)


  #Modèle AdapSCAD:
var_ascad <- df_beta %>% select(1,11) %>% filter(AdapSCAD != 0) %>% select(1)

mod_ascad <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Container.Index + WTI + CPI + Indice_Kilian + Three_Component_Index + News_Based_Policy_Uncert_Index + world.trade..volume. + Business_Tendency_Surveys + Capacity_Utilization + Consumer_Sentiment + GISS_Temperature + GEPU_ppp + Industrial_Capacity + OECD_6NME_Industrial_Production + Spread + Trade_Weighted_USD, data = df_fret_clean)



 #Modèle MS_AdapEN

var_msaen <- df_beta %>% select(1,12) %>% filter(MS_AdapEN != 0) %>% select(1)

mod_msaen <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Indice_Kilian + Three_Component_Index + Consumer_Sentiment + GEPU_ppp + Trade_Weighted_USD, data = df_fret_clean)



 #Modèle MS_AdapSCAD   (seulement News_Based_Policy_Uncert_Index pas selectionné comparer à AdapSCAD)

var_msascad <- df_beta %>% select(1,13) %>% filter(MS_AdapSCAD != 0) %>% select(1)

mod_msascad <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Container.Index + WTI + CPI + Indice_Kilian + Three_Component_Index + world.trade..volume. + Business_Tendency_Surveys + Capacity_Utilization + Consumer_Sentiment + GISS_Temperature + GEPU_ppp + Industrial_Capacity + OECD_6NME_Industrial_Production + Spread + Trade_Weighted_USD, data = df_fret_clean)
 
 #Modèle BoLASSO 

mod_bolasso <- lm(Baltic.Dry.Index ~ lag(Baltic.Dry.Index) + Gold_Princing_Index + News_Based_Policy_Uncert_Index + Capacity_Utilization + GISS_Temperature + Industrial_Capacity + OECD_6NME_Industrial_Production + US_Ending_Stocks_Crude_Oil + VIX, data = df_fret_clean)

Tableau récapitulatif

Resultat des estimations
	GETS	LASSO	Elastic-Net	WF_LASSO	Bridge	SCAD	AdapLASSO	AdapEN	AdapSCAD	MS_AdapEN	MS_AdapSCAD	BoLASSO
(Intercept)	-7.708	-0.611	15.597	7.633	5.127	66.453	6.414	15.597	58.792	15.597	54.866	-38.228
	(62.444)	(44.666)	(41.060)	(40.869)	(43.255)	(50.318)	(43.219)	(41.060)	(67.556)	(41.060)	(67.121)	(68.245)
lag(Baltic.Dry.Index)	-0.099	-0.073	-0.124	-0.109	-0.131	-0.083	-0.110	-0.124	-0.132*	-0.124	-0.128	0.264***
	(0.076)	(0.081)	(0.075)	(0.075)	(0.082)	(0.076)	(0.079)	(0.075)	(0.079)	(0.075)	(0.078)	(0.081)
Container.Index	122.616								80.178		83.937
	(74.662)								(83.546)		(83.145)
CPI	-172.251**					-147.472*			-267.077***		-263.296***
	(81.096)					(79.705)			(94.194)		(93.789)
Indice_Kilian	19.279***	18.356***	18.898***	18.732***	18.001***	19.138***	18.184***	18.898***	19.045***	18.898***	19.058***
	(2.161)	(2.365)	(2.192)	(2.177)	(2.296)	(2.173)	(2.287)	(2.192)	(2.300)	(2.192)	(2.295)
Consumer_Sentiment	-29.338***		-28.490**	-26.961**		-28.597**		-28.490**	-28.505**	-28.490**	-28.734**
	(11.172)		(11.428)	(11.368)		(11.234)		(11.428)	(11.649)	(11.428)	(11.616)
GEPU_ppp	-3.176*		-3.399*	-2.834		-2.938		-3.399*	-3.616*	-3.399*	-3.458*
	(1.909)		(1.909)	(1.935)		(1.919)		(1.909)	(1.999)	(1.909)	(1.980)
Industrial_Capacity	549.267*								523.194		533.214	337.095
	(326.734)								(336.401)		(335.250)	(411.627)
Trade_Weighted_USD	-112.289***		-85.546***	-91.686***	-89.525**	-113.327***	-106.087***	-85.546***	-77.540**	-85.546***	-78.814**
	(32.898)		(31.309)	(31.398)	(36.654)	(32.823)	(32.246)	(31.309)	(36.464)	(31.309)	(36.325)
Three_Component_Index	-3.629*		-3.684*	-3.566*		-3.803**		-3.684*	-8.608	-3.684*	-3.388*
	(1.902)		(1.947)	(1.936)		(1.914)		(1.947)	(8.353)	(1.947)	(1.944)
US_Ending_Stocks_Crude_Oil	0.006*			0.005		0.005						0.006
	(0.004)			(0.004)		(0.004)						(0.005)
Petroleum_and_other_liquids_stocks_US				0.003
				(0.003)
WTI					8.114				17.278*		17.028
					(8.528)				(10.340)		(10.309)
News_Based_Policy_Uncert_Index									3.237			-1.907
									(5.037)			(1.368)
world.trade..volume.									-35.928		-36.064
									(40.915)		(40.818)
Business_Tendency_Surveys									114.820		107.029
									(234.118)		(233.254)
Capacity_Utilization									-5.479		-5.334	204.475
									(121.667)		(121.381)	(147.459)
GISS_Temperature									444.159		470.794	482.459
									(375.717)		(372.547)	(481.382)
OECD_6NME_Industrial_Production									3.436		4.534	-61.583
									(84.706)		(84.490)	(94.497)
Spread									-218.195		-215.996
									(246.281)		(245.678)
Gold_Princing_Index												1.841
												(1.214)
VIX												-23.913
												(15.137)
Num.Obs.	142	142	142	142	142	142	142	142	142	142	142	142
R2	0.513	0.354	0.472	0.490	0.405	0.498	0.401	0.472	0.526	0.472	0.524	0.155
R2 Adj.	0.476	0.344	0.449	0.459	0.387	0.468	0.388	0.449	0.461	0.449	0.463	0.098
AIC	2166.4	2190.5	2169.7	2168.9	2182.9	2166.6	2181.8	2169.7	2176.6	2169.7	2175.1	2242.6
BIC	2201.9	2202.4	2193.4	2198.5	2200.6	2196.2	2196.6	2193.4	2232.8	2193.4	2228.3	2275.1
Log.Lik.	-1071.192	-1091.275	-1076.872	-1074.474	-1085.446	-1073.296	-1085.914	-1076.872	-1069.315	-1076.872	-1069.551	-1110.305
F	13.798	38.045	20.147	15.969	23.284	16.514	30.764	20.147	8.084	20.147	8.604	2.693
* p < 0.1, p < 0.05, * p < 0.01

On constate ainsi, à partir du tableau récapitulatif des modèles estimés, que la régression linéaire la plus optimale semble être celle où ont été incluses les variables retenues par la méthode GETS. La qualité d’ajustement est effectivement la plus élevée avec \(R^2\)=47.6%, l’AIC le plus faible (2166.4). Pour rappel ce modèle considérait les variables “Container.Index”, “CPI”, “Indice_Kilian”, “Consumer_Sentiment”, “GEPU_ppp”, “Industrial_Capacity”, “Trade_Weighted_USD”, “Three_Component_Index”, “US_Ending_Stocks_Crude_Oil” comme importantes dans l’explication du BDI. Par rapport aux autres méthodes de sélection de variables, la régression GETS avait l’avantage d’en considérer un plus grand nombre tout en restant parcimonieuse notamment en comparaison avec le LASSO qui n’en gardait qu’une seule ou bien l’adaptive SCAD qui en sélectionnait 16. Cependant nous devons vérifier les hypothèses sur les résidus avant de pouvoir conclure et interprêter les résultats.

On peut aussi faire le constat que l’ajout de la variable à expliquer laguée en régresseur n’améliore pas les estimations puisque la variable n’est significative dans aucune des estimations exceptée l’AdaptSCAD au seuil de risque de 10% et le BoLASSO à 1%. Par ailleurs, le modèle le plus pertinent qui arrive en 2ème ‘position’ par les critères d’informations, de qualité et de vraisemblance, est le modèles SCAD.

L’objectif de ce dossier était de tester différentes méthodes de sélection de variables étudiées en cours, mais aussi d’en chercher une en plus dont il fallait expliquer le fonctionnement avant de l’appliquer. C’est ce que nous avons fait pour la méthode du LASSO bootstrapé, mais nous avons montré comment son application était compliquée du fait des packages non valables pour notre version R studio, ni sur R, ou alors qui étaient prévus pour des natures de Y différentes de la nôtre. Aussi nous avons utilisé le package “mhtbooot” qui se rapprochait du package “mht” qui contenait la fonction bolasso() pour l’application de notre méthode. Cependant, au vu des variables sélectionnées qui ne correspondent à aucune autre sélection et au vu de la qualité d’ajustement de la régression linéaire comportant les 8 variables sélectionnées par BoLASSO, nous comprenons que le package ne doit pas implémenter exactement la fonction que nous recherchions. Par conséquent nous laissons cette méthode de côté et n’effectuerons même pas l’analyse des résidus sur le MCO correspondant. En revanche théoriquement cette méthode est pertinente et il serait intéressant de creuser son application en utilisant peut-être un autre logiciel d’analyse.

Par ailleurs, nous voyons que les résultats sont identiques pour les modèles composés des variables retenues par les méthodes de sélection “Elastic-Net”, “Adaptive EN” et “Multi-Step Adaptive EN”, ce qui est logique puisque comme nous pouvons le constater, ces 3 méthodes ont sélectionné les mêmes régresseurs du BDI.

Enfin, nous passons à l’analyse des résidus de chacun des modèles pour la validation de ceux-ci.

Analyse des résidus

Chaque modèle estimé par la méthode des moindres carrés ordinaires doit répondre à des critères spécifiques que nous décrivons ci-dessous, avant que le modèle puisse être validé et que des conclusions puissent être tirées.

Le modèle doit être bien spécifié, c’est-à-dire qu’il existe une relation linéaire entre les variables, entre les paramètres et qu’il n’y a pas d’omission ou de redondance des variables. Le test Reset Ramsey nous permettra de vérifier si notre modèle est bien spécifié ou non. Il ne doit pas non plus y avoir de multicollinéarité parfaite entre les variables exogènes, en d’autres termes, le coefficient de corrélation simple entre deux variables explicatives doit être différent de l’unité. Nous le vérifierons avec le Variance Inflation Factor, VIF.

Les dernières hypothèses concernent les perturbations aléatoires. Premièrement, les résidus doivent suivre une distribution normale. Ce qui peut être vérifié avec le test de Jarque-Bera mais ne doit être vérifié que lorsque l’échantillon comporte moins de 30 observations. Dans notre cas, nous avons 142 observations, il ne sera donc pas nécessaire de vérifier cette hypothèse. Deuxièmement, l’homoscédasticité des résidus doit être contrôlée c’est-à-dire que la variance du terme d’erreur doit rester constante pendant toute la période de l’échantillon. Plusieurs tests existent pour cette hypothèse, mais nous utiliserons ici le test de Breusch-Pagan. Si toutes ces conditions ne sont pas vérifiées, le modèle ne sera pas valide.

	gets_beta	LASSO	Elastic_Net	WF_LASSO	Bridge	SCAD	AdapLASSO	AdapSCAD	MS_AdapSCAD
Reset_Test	0.003687672	0.697745732	0.132101945	0.030430639	0.978847841	0.006174374	0.937414194	0.002132790	0.002876209
Breush_Pagan_Test	0.19890757	0.08227370	0.63398549	0.89419397	0.19488943	0.77290191	0.21269083	0.01597885	0.01113370
Vif_test	1.548846	1.423771	1.456205	1.486587	1.565204	1.529550	1.453179	26.544930	2.569443
Conclusion	X	Ok	Ok	X	Ok	X	Ok	X	X

Les résultats des tests sur la validité des conditions d’un MCO apparaissent dans la table ci-dessus. Pour les 3 modèles avec les mêmes variables retenues par les méthodes de sélection de variables (“Elastic-Net”, “Adaptive EN” et “Multi-Step Adaptive EN”) nous ne procédons évidemment aux tests de vérification que pour un seul d’entre eux que nous laissons sous le nom de “Elastic-Net”.

Concernant la forme fonctionnelle du modèle et comme évoqué précédemment, c’est le test de Ramsey qui nous informe de sa bonne spécification, il est basé sur les hypothèses suivantes :

\(H_0\) : la forme fonctionnelle du modèle est linéaire
\(H_1\) : la forme fonctionnelle du modèle n’est pas linéaire

On voit ainsi que 4 modèles sur 9 valident cette hypothèse puisque la p-value associée au test Reset est supérieure au seuil de risque de 5%. Le modèle GETS qui paraissait le mieux d’un point de vue de la qualité d’ajustement ne valide pourtant pas cette hypothèse et a donc une forme fonctionnelle non linéaire.

Le test de Breusch-Pagan quant à lui, qui concerne l’homoscédasticité des résidus suit la règle de décision suivante :

\(H_0\) : les résidus du modèle sont homoscédastiques
\(H_1\) : les résidus du modèle sont hétéroscédastiques

On constate cette fois que seuls les modèles avec les variables sélectionnées par le Multi-Step Adaptive SCAD et l’adaptive SCAD ont des résidus hétéroscédastiques. Tous les autres valident l’hypothèse fondamentale à l’acceptation d’un modèle de régression linéaire.

Enfin, la vérification de la multicolinéarité au sein des variables explicatives nous est donnée par le critère du “Variance Inflation Factor” (VIF) qui, bien que ce ne soit pas un test, s’applique de la manière suivante :

si le VIF est supérieur à 10 alors il y a un problème de multicolinéarité

Nous avons mis pour cette condition le VIF le plus élevé parmi les variables explicatives ; nous voyons que même la valeur la plus grande est loin d’atteindre le seuil critique de 10 où nous sommes contraints de retirer la variable sauf pour le modèle adaptive SCAD. Le vif le plus élevé est en effet de 26.54 correspondant à la variable “News_Based_Policy_Uncert_Index”.

Par conséquent, bien que le meilleur modèle en termes de qualité soit celui où ont été incluses les variables issues de la méthode de sélection GETS, nous voyons que le régression ne valide pas les hypothèses fondamentales de validité et ne peut donc pas être retenue. Ainsi, si l’on regarde le meilleur modèle après celui-ci (i.e. celui qui est le plus pertinent ; \(R^2\), AIC, Log.Llh et qui valide toutes les hypothèses) on voit qu’il s’agit de celui où ont été incluses les variables sélectionnées par les méthodes “Elastic-Net”, “Adaptive EN” et “Multi-Step Adaptive EN”. Leur qualité d’ajustement est de 44.9% et leur critère Akaike de 2169.7.

Nous pouvons donc interpréter les résultats de la régression linéaire prenant en compte les variables retenues par ces 3 méthodes. Attention cependant aux valeurs des coefficients car nous avons stationnarisé toutes nos variables, donc les coefficients apparents sont ceux des rentabilités. On voit que le modèle est composé de 6 variables explicatives dont 3 sont significatives au seuil de risque 5% et 2 le sont uniquement au seuil de 10%.

Ainsi, lorsque les rentabilités de l’indice Kilian augmentent d’une unité, les rentabilités du BDI vont augmenter elles aussi de 18.9 points. La variable qui a l’effet le plus grand sur Y est celle du commerce américain (“Trade_Weighted_USD”) bien que, n’étant pas un indice boursier, nous ne pouvions interpréter sa valeur comme une rentabilité. Par ailleurs, il apparaît que lorsque les rentabilités du “Three Component Index” augmentent d’une unité celles de l’indice maritime diminuent de 3.68 points, au seuil de risque 10%.

Ainsi, l’application des méthodes de sélection de variables a été concluante sur notre base puisqu’elles ont permis la construction de modèles parcimonieux et explicites. Nous avons alors pu expliquer 44.9% des variations des rentabilités du Baltic Dry Index sur la période février 2017 - décembre 2018.

5- Annexes

Annexe n°1 : Graphiques des séries en niveau

Annexe n°2 : Graphiques des séries transformées

Annexe n°3 : Distribution des variables brutes et transformées

Annexe n°4 : Pourcentage de la variance expliquée par les axes, ACP

Prévisions du Baltic Dry Index et méthodes de sélection de variables

Diane THIERRY et Teodoro MOUNIER TEBAS

2020-12-05

Sommaire

1. Introduction

2. Analyses exploratoire

2.1 Importation

2.2 Analyse de la stationnarité

2.3 Analyse des outliers

2.4 Graphiques

2.5 Statistiques descriptives

2.6 Méthodes de classification

2.7 Corrélation

3. Sélection des variables

3.1 Approche économétrique : GETS

Théorie

Application

Fonction ‘isat’

3.2 Régressions pénalisées :

3.2.1 Ridge regression

Théorie

Application

3.2.2 Lasso regression

Théorie

Application

3.2.3 Weighted fusion Lasso regressions

Théorie

Application

3.2.4 Elastic-Net regression

Théorie

Application

3.2.5 Bridge regression

Théorie

Application

3.2.6 SCAD

Théorie

Application

Complément

3.2.7 Adaptive Lasso regression

Théorie

Application

3.2.8 Adaptive Elastic-Net

Théorie

Application

3.2.9 Adaptive SCAD-Net

3.2.10 Multi-Step Adaptive Elastic-Net

Théorie

Application

3.2.11 Multi-Step Adaptive SCAD-Net

3.3 Autre approche :

Théorie

Application

3.4 Comparaison

4. Modèle de régression

Code des modèles

Tableau récapitulatif

Analyse des résidus

5- Annexes